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Non-separable bidimensional 
wavelet bases 

Albert Cohen and Ingrid Daubechies 



Abstract, We build orthonormal and biorthoepa] wavelet bases of 
X 2 (R 2 ) with dilation matrices of determinant 2. As for the one dimen- 
sional case, oufconstruction uses a scaling function which solves a two- 
scale difference equation associated to a Flit filter. Our wavelets are 
generated from a single compactly supported mother function. How- 
ever, the regularity of these functions cannot be derived by the same 
approach as in the one dimensional case. We review existing techniques 
to evaluate the regularity of wavelets, and we introduce new methods 
which allow to estimate the smoothness of non-separable wavelets and 
scaling functions in the most general situations. We illustrate these 
with several examples. 

I. Introduction. 

. In the most general sense, wavelet bases are discrete families of 
functions obtained by dilations and translations of a finite number of 
well chosen mother functions. The most well known are certainly dyadic 
orthonormal bases of L 2 (I), of the type 

(1.1) i>{{x) = 2-^(2^-1), f,KZ. ' 

These constructions have found many interesting applications^ both in 
mathematics because they form Riesz bases for many functional spaces' 
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and in signal processing because wavelet expansions are more appropri- 
ate than Fourier series to represent the abrupt changes in non-stationary 
signals. 

Several examples have been given by Meyer [Mel], Lemarie [Le] 
and Daubechies [Daul], generalizing the classic Haar basis in which 
the mother wavelet </> = X[0)1/2] - X[i/ 2) i] suffers from a lack of 
regularity.since it is not even continuous. All afe based on the concept 
of multiscale analysis, it. a ladder of closed subspaces {Vj} m which 
approximates I 2 (t), 

(1-2) {0}->...y 1 c7 0 CF_ 1 ...->I 2 (R), 

(note that in some papers and in Meyer's book, the converse convention 
is used, it. Vj C %i) and satisfies the following properties, 

(1.4) there exists a function <p{x) in V 0 such that the 
set {tp(x - k)}kzz is an orthonormal basis for 7 0 . 

Since 7 0 C V- h the scaling function ip(x) has to be the solution of a 
two scale difference equation, 

(1.5) ^) = 2^ Cn ^-n). 

The associated wavelet is then derived from the scaling function by the 



M *M = 2£(-l)'? Hiy (fr-n). 

na 

In the standard interpretation of a multiresolution analysis, the 
projections of a function / on the spaces Vj are viewed as successive 
approximations to /, with liner and finer resolution as j decreases. The 
wavelets can then be used to express the additional details needed to 
go from one resolution to the next finer level, since the {ij)(x - k)} k& 
constitute an orthonormal basis for W 0 , the orthogonal complement of 
Vo in V-l The whole set {^(j) } j ka forms then an orthonormal basis 
ofL 2 (E). 

: We are here interested in similar constructions adapted to functions 
or signals of more than one variable. 
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„ The most commonly used method to build a multiresolution anal- 
ysis and wavelet bases in L\l n ) is the tensor product of a multires- 
olution analyses of'tffR). In £ 2 (1 2 ) it leads to a ladder of spaces 
Vj = Vj ® Vj C Vj-i generated by the families, 

(1.7) ♦L(*,y) = r^i-t) y (r^-!), Mez: 

Three wavelets are then necessary to construct the orthogonal comple- 
ment of Vo in V_i, namely, 



(1.8) 
(1.9) 
(1.10) 



%2/)=#My), 



Actually, the theory of multiresolution analysis, as it was intro- 
duced by S. Mallat and Y. Meyer (see [Mai] and [Mel]) was first mo- 
tivated by the possibility of building these separable wavelets for the 
analysis of digital picture. 

It is clear,'however, that this choice is restrictive and that it gives a 
particular importance to the x and y directions, since $ tt and $j match 
respectively the horizontal and vertical details. 

A more general way of extending multiresolution analysis to n di- 
mensions consists in replacing the axioma (1.3) and (1.4) by 

(1.11) ; '/(,) 6 V, *=►./(&) £ V H 

(1.12) There exists a function j in V 0 such that the set 
{<l>(x - fc)} ieZ n is an orthonormal basis for V 0 , 

g- i where D is a n x n dilation matrix. 

All the singular values Ai , . . . , A n of D must satisfy 



(1.13) 



W>i, 



to ensure that the approximation gets finer in every direction as j goes 
to -oo. Furthermore, we require D to have integer entries. This con- 
dition means that the action of D on the translation grid Z" leads to a 
sublattice Y C Z n . 

The number of basic wavelets required to characterize the orthog- 
onal complement of V 0 in V.\ is in that case trivially given by the 
following heuristic argument. This complement should be generated « 
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by the action of Z" on the basic wavelets, in the same way that \ is 
generated by the action of Z n on A whereas Y A is generated by the 
action ofD^Z". Consequently, each of the generating functions can be 
associated with an elementary coset of d~ l l n jl n ~ Z n /DZ n except 
one which corresponds to the scaling function (see figure 1). There- 
fore, d = |detD| - 1 different wavelets are needed. Note that it is 
not strictly necessary that the entries of D be integer to build wavelet 
bases using D as the elementary dilation, 



-U 



(-1,0) 



(M) 



Figure 1 

Z 2 and DZ 2 in the case where D - 



2 -1 
1 2 



The scaling function and the four basic wavelets 
are indexed by an element of Z 2 /DZ 2 . 
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However, the condition seems to be necessary for the existence of a mul- 
' tiresolution analysis based on a single, real valued, compactly supported 
scaling function. 

In this work we shall indeed focus on real valued, compactly sup- 
ported scaling functions and wavelets. They have the advantage that 
the sequence {c n } n€ z introduced in the two scale difference equation 
(1.5) is real and finite. These coefficients play an important part in 
the numerical applications because they are used directly in the Fast 
Wavelet Transform algorithm as decomposition and reconstruction fil- 
ters. They constitute in that case an FIR (finite impulse response) filter 
which can be implemented very easily. Furthermore, this finite set of 
coefficients contains all the information about the multiresolution anal- 
ysis since the functions <p and j can be constructed as solutions of (1.5) 
and (1.6). Our starting point to build wavelet bases will thus be afinite 
set of coefficients and the associate two-scale difference equation, rather 
than the approximation spaces Vj themselves. 

The main difficulty in this approach is the design of the FIR filter 
{ c n}n=o,.,w in such a way that <p and j> are smooth and have orthonor- 
mal translates. 

In the one dimensional case, it is shown in [Daul] that orthonor- 
mal wavelets can be constructed by choosing a filter which corresponds 
to a particular case of exact reconstruction subband coding schemes, 
and which can be made arbitrarily regular by increasing the number 
of taps in a proper way. Several contributions have followed, giving 
supplementary information on the type of filter which has to be used 
(see [Me2],'[DL], [Col], [Dau2], [Co2J, [Dau3]). 

In the present bidimensional case, the design of filters associated 
to "nice wavelet bases" turns out to be more difficult because some of 
the one-dimensional techniques do not generalize trivially (or do not 
generalize at all!) to higher dimensions and new methods have to be 
introduced, This article concentrates on the situation where D is a 2 x 2 
matrix with | detD| = 2. 

We deliberately restrict ourselves to this set of matrices for two 
reasons: 

• These dilations have already been considered by electrical engi- 
neers and seem to have interesting applications in signal analysis 
and image processing. For example, since only one basic wavelet 
is required, one may hope for a more isotropic analysis than with 
the separable construction. Subband coding schemes with -deci- 
mation on the quincunx sublattice have been studied in the works . 
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of J. C. Feauveau [Fea] and M. Vetterli and J, Kovacevic [KV], 
Our work is complementary to their signal processing approach 
since we investigate here the mathematical properties, such as the 
Holder regularity of the wavelet bases associated to these schemes. 
This regularity is important when one asks that the reconstruc- 
tion of the signal from the coarse scales has a smooth aspect (see 
Section II.2). 

• These dilations are simple and our study will be reduced to the 
case of two basic matrices. However, the difficulties which appear 
in the evaluation of the regularity of the corresponding wavelets 
are common to all the non-separable constructions, and the tech- 
niques that we develop to solve this problem can be used for other 
types of dilations. We believe that the set of integer matrices with 
| det X) | = 2 constitutes an interesting "laboratory case" in the 



In the next section of this paper, we shall give an overview of 
different techniques which can be used in the construction of one di- 
mensional compactly supported wavelets. Some new tools will be intro- 
duced specifically to be generalized and used in the multidimensional 
situation. 

The third section examines the possible subband coding schemes 
with decimation on the quincunx sublattice and their general relations 
with non-separable wavelet bases. 

In the fourth section, orthonormal bases of wavelets are constructed 
from such coding schemes. We show that for the same filters, different 
bases with widely differing regularity can be obtained, depending on the 
choice of the dilation matrix. Finally, we use a biorthogonal approach, 
in Section V, to construct more symmetrical wavelet bases correspond- 
ing to linear phase filters and allowing a more isotropic analysis. We 
show that arbitrarily high regularity can be attained and we give some 



II. The construction of compactly supported wavelets in one 
dimension: A complete toolbox, 

The purpose of this section is to review, in the one dimensional 
case, many different techniques that can be used to build regular wave- 
lets from subband coding schemes, theoretically and numerically. Some 
of these techniques, like the Littlewood-Paley estimation of smoothness, 
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are not frequently used In the one dimensional case,. but they turn out 
to be very useful for the non-separable bidimensional wavelets For 

wtB, can also consult [M M !Ma11, M 



Wavelet bases and subband coding schemes. 
Il.l.a. The orthonormal case. 

Let {Vj} i& be a multiresolution analysis of X 2 (R). We can use 
the discrete Fourier transform of the finite sequence {c n }^ w , U the 



(2-1) m 0 ( W ) = £ Cn e-- = £ Cn e--, 

to rewrite the two scale difference equation (1.5) that characterizes ip(x). 
We suppose that the c n are real. Taking the Fourier transform of (1.5) 
and (1.6) we obtain 

(2.2) £(2w) = m o (w}0(w) 



(2.3) ^(2o;) = e- iw mo(a; + 7r)^) = m 1 (^( w ). 

Two fundamental properties of m 0 (u>) can be derived from the mul- 
tiresolution analysis properties 

• Since {^(2; - k)} k& is an orthonormal basis of Pq, the Fourier 
transform <j>(w) satisfies a Poisson identity 



(2.4) £|^ + 2mt)| 2 = l. 

Combined with (2.2) this leads to 

(2.5) |moH 2 •+ |m 0 (o; + 7r)| 2 = l 
which may also be written as 

(2.6) 2 £ c n c ni2k = f M (= 1 if h = 0, 0 otherwise) . 

nez j 
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• The denseness of {Vj}^ in I 2 (R) is equivalent to 

#(0) = j<p{x)ix = 1, 

(see [Mel], [Mai] or [Col]). 
Consequently, we have 

(2-7) m 0 (0) = 1 and m^) = 0, 

which may also be written as 

^2 jV 2 

(2-8) ' ^c n =:l and £(-l)« Cn = 0 . 

n=Wi n=Ni 

The subband coding scheme associated to our multiresolution analysis 
appears clearly in the Fast Wavelet Transform Algorithm of S. Maljat 
[Ma2]. Let us recall how it works. The initial data are considered as 
the approximation of a continuous function at the scale j - 0, 



(2.9) 



This allows the computation of the approximations and the details at 
coarser scales, i.e. 

(2.10) S>=Vl 2 (/yj and D[ = T» (/.^ j>0. 



(The coefficients are normalized in such way that if / = 1 locally, then 
S\ = 1 in that area). The sequence {S[} ka (respectively {Dj}^ 
is then derived from {S\~% a by a convolution with the filter mJ(J 
(respectively, m^w)) followed by a decimation of one sample out of two 
to keep the same total amount of information, it. 



The algorithm then iterates on {S[}keh Conversely, the sequence 
{S[ }fc € z can be recovered by applying the same filters m 0 (w) and 
mi(w) on {5j}jbgz and {Djjjtgz after inserting a zero between every 
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pair of consecutive samples, and summing the two 



i, i.e. 



All these operations, decomposition - decimation - interpolation - re- 
construction, constitute a complete subband coding scheme as shown 
on figure 2. The property of exact reconstruction can now be derived in 
two ways. It is a natural consequence of the multiresolution approach, 
since Vj = Vfa e W m but it can also be viewed as a consequence 
of formula (2.5) for the filter m 0 . This type of filter pair (m 0) mA is 
known as a pair of "conjugate quadrature filters" (CQF); they were 
first discovered by Smith and Barnwell in 1983, c/. [SBl], The design 
of FIR pairs, with real coefficients and perfect reconstruction, has been 
generalized in [Daul], It also appears in [ASH], [SB2], [Vel].' 




reconstructed 



Figure 2 

Subband coding scheme corresponding to the FWT algorithm. 
The sign 2 j stands for "decimation of one sample out of two" and 2 j for 
the insertion of zeros at the intermediate values. 

• Since m 0 (w.) is regular (it is a trigonometric polynomial) and since 
m 0 (0) = 1, we can iterate (2.2) to obtain 



(2.11) 



+00 



Given a conjugate quadrature filter m 0 (w) (i.e. a trigonometric polyno- 
mial satisfying (2.5) and (2.7)), it is thus possible to define the scaling 
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1, either as a solution of the two s 
or explicitly with the above infinite product. However, this does not 
always lead to a multiresolution analysis: the function ip(x) = |x[o 3] 
generated by the CQF m 0 (w) = (1 + e 3iw )/2, for example, does not 
satisfy the orthonormality of the translates. Orthonormality of the 
ip(x - k) turns out to be equivalent to the L 2 convergence of the trun- 
cated products 0„(w) = nj sl m^u)^^) to ftp) (be- 
cause {ip n (x-k)} k& is an orthonormal set as soon as (2.5) is satisfied). 

More precisely, the following result characterizes the subclass of 
CQF filters leading to a multiresolution analysis and orthonormal basis 
of wavelets. 

Theorem 2.1. Let m 0 (w) he a Conjugate Quadrature Filter. Then, 
the infiniie yroiud (2.11) his to a multiresolution analysis ij ani only 
if there exist a compact set Kct such that, 

i) K tonkins a neighhourhooi of the origin, 

ii) \K\ - It mi for all u m.[-7r, t], there exist n 6 Z such that 

iii) for all n > 0, m 0 (2" n w) dots not vanish on K. 

The set /{ is said to be "congruent to [-tt, jt] modulo 27r" (figure 
3). The proof of this result can be found in [Col]. It exploits the 
continuity of mo, the compactness of K and mo(0) = 1 to show that 
(iii) is equivalent to p(w) > c > 0 on K. This is then sufficient to derive 
the L 2 convergence of the % by Lebesgue's Theorem. We shall use a 
multidimensional generalization of Theorem 2.1 in the fourth section. 
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Figure 3 

Example of compact set congruent to [-7T, t] modulo 2f. 
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ILl.b, The biorthogonal case. 

The conjugate quadrature filters are a very particular case of sub- 
band coding scheme with perfect reconstruction, because identical fil- 
ters (up to a complex conjugation) are used for both the decomposition 
and the reconstruction stages. If we do not impose this restriction, 
•then the scheme uses four different filters: m 0 (w) and m^w) for the 
decomposition, m 0 (w) and m^w) for the reconstruction. Perfect re- 
construction for any discrete signal is then ensured if, 



(2.12) 



m 0 (w) m 0 (w) + mi(w) roi(w) = 1 « 

( m 0 (« + tt) m 0 (w) 4 mi(w + jt) m$)-. = 0 . 

m 0 (w) and mi(w) may thus be regarded as the solutions of a linear 
system.. However, to avoid the infinite impulse response solutions, we 
' shall force the determinant of this system to be ae'* w , a f 0, k 6 Z. For 
sake of convenience we take a = -1 and k ~ 1 (a change of these values 
would only mean a shift and a scalar multiplication on the impulse 
response of our filters). This leads to 

(2.13) m 0 (w) m 0 («) + rao(w + ir) m 0 (w + f ) = 1 , 



(2.14) mi(«) = e^m 0 (w+ir) > m : {u) = f iu m 0 (o; + 7r) . 

The formulas (2.13) and (2.14) are thus the most general setting for 
finite impulse response subband coders with exact reconstruction (in 

■ the two channel' case). Thefunctions m 0 (w) and m 0 («) are called "dual 
filters". It is clear that the special case m 0 (w) = m 0 (w) corresponds 
to the conjugate quadrature filters of Il.l.a. However, dual filters are 
easier to design than CQF's. For example, if m 0 ' is fixed; m 0 can be 
found as the solution of a Bezout problem which is equivalent to a linear 
system. The coefficients of; these filters can be very simple numerically 
(in particular they can have finite binary expansion which is very useful 
for practical implementation), furthermore they can be chosen symmet- 
rical ("linear phase filter"), a property which is impossible to satisfy in 
the CQF case. ; j . 

We can mimic, in this more general framework, the construction of 

■ orthonormal wavelets from CQF. Assuming that m 0 (0) = m 0 (0) = l 
andm'o(Tr) = m 0 (7r) = 0, we define 



(2.15) 



W = ]] mo(2~ A a;), 
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(2.16) 



i{2w) = roi(w)$w) , 



(2.17) fa) = n*o(2"M, 

(2.18) j(2u) = mi(w)^(«) . 
In [CDF], the following theorem was proved, 
Theorem 2.2. 

,J / Ww) = nLim 0 (2- fc w)x[-2" ff ,2^]M mi Ww) = K=i 
m 0 (2~*w) X[-2"7r,2 n 7r](w) wwerye in £ 2 (R) rcapeciwefy to #(w) 
ani ij>(w) } then the following duality relations tire satisfied ' 

(2.20) (H^^kv 
mi for all f in L 2 (M) one has the wipe decoropoafc 



(2.21) 



j=-j ket 



(in the I) sense). 

• If ip ani (p satisfy \<j>(w)\ + |p(w)| < C (1 + |^|)~ 1/2_fi: /or some 
e > 0, to the families {^J^gz flwJ W{W fi?, e jraroej o/ 
L 2 (E). 

• WTieri ieae too properties holi, then ^JJ^gi ore fiior%- 
ono/ (or kal)Rie3z oases of L 2 (I). 

Many examples of these systems can be found in [CDF] and a 
sharper analysis of the frame conditions is developed in [CD]. We now 
recall a practical way of constructing ip and $ numerically from a given 
subband coding scheme. 



II.2. The cascade algorithm. 

In the last section we saw that the scaling function <p(x) could be 
approximated, at least in L 2 (E), by a sequence of band limited functions 
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M„>o defined by 



(2.22) ' ^) = Umo(2^) X [-2« ffl ^]M. 
i=i 

These functions are characterized by their sampled values "at the points 
2-Jfe(i€Z),i.e. 



(2.23) 



s n k = 9n(2~ n k). 



This sequence can also be considered as the impulse response of the 
transfer function 



(2.24) S n (u) = 2 n ][m^). 

S n (u) can be obtained recursively by the formula 

(2.25) 5 B+1 ( W ) = 2mo(«)S»(2w). 

In the time domain, (2.25) becomes an interpolation scheme; the se- 
quence aj is dilated by insertion of zeros (S n (u) -> S n [2u)) before 
being filtered (multiplication by 2m 0 (w)). We have thus, 

(2.26) s^^lc^sl IT" 

This iterative process, which computes the {sj}^ sequences from w 
initial Dirac sequence ft,* is called the "cascade algorithm". We illus- 
trate it on figure 4 (our sequences are represented by piecewise constant 
functions). 

Note that it identifies exactly with the reconstruction stage in the 
FWT algorithm described in Il.l.a. The scaling function is thus 
approached by the reconstructed signal from a single approximation 
coefficient at a coarse scale. Similarly, the wavelet will be obtained bj 
starting the reconstruction from a detail coefficient at a coarse scale 
(and thus applying mi (id) at the first step of the cascade). 

This explains why subband coding schemes associated with regu- 
lar wavelets are particularly interesting: the smoothness of the wavelel 
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determines the appearance of the coarse scale components of the recon- 
structed signal. A smooth appearance is important for many applica- 
tions such as compression where a big part of the finer scale information 
is thrown away. 

In the biorthogonal case, the analysis and the synthesis wavelets [ij) 
and fj need not have the same regularity. As just discussed, smoothness 

«) h 0 1 



_ ,i 

-1/2 0 1/2 




-! -1/2 0 1/2 1 




Figure 4 

The cascade algorithm (from [Daul]). 
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is important for the reconstructing function; the analyzing function 
needs only to be sufficiently regular to ensure that the wavelet bases 
are unconditional, so that the FWT algorithm is stable. Note that 
an important property on the analyzing wavelet is cancellation, i.e. 
vanishing moments, ensuring small high scale coefficients for smooth 
regions in the function or signal to be analyzed. 
|| Let us finally mention that this type of "refinement method" is 

well known in approximation theory as "stationary subdivision" (e.^. 
;f , [CDM], [DyL]). Most of these papers are motivated by interpolation 
problems, where smooth curves or surfaces need to be constructed, con- 
necting (or close to) given sparse data points. Consequently, they are 
mainly concerned with what we call the reconstruction stage and they 
do not study the existence of an associated subband coding scheme. 
This also means that they do not care about an easy way of encoding 
; I or representing, the extra "detail information" (<— > Wj) that xan be 
iff added in going from one refinement level to the next one (Vj ->7;-i)- 
j |; On the other hand, the subband coding literature seldom mentions the 
importance of the smoothness appearing in the cascade of the recon- 
struction from the low scales. Orthonormal and biorthogonal wavelet 
bases lead to an elegant combination of these two approaches. 
:|. - , . We now present several different methods to estimate the regular- 
I ity of the wavelets associated to a given subband coding scheme. We 
shall concentrate on the regularity of the scaling function which deter- 
| J . mines the regularity of the wavelet itself because t/>(z) is a finite linear 
combination of translates of <p(2x). Whatever the method used, if a 
\ |; global regularity of order r is achieved, then the cascade algorithm also 
\\ converges uniformly up to this order (see [Daul], [DL], [Co2]). 

; : f IL3. Regularity: the spectral approach. 

II.3.a. A Fourier estimation of the Holder exponent. 

■ 

Let us denote by C a the Holder space defined as follows. For 
I . a = n-H0,/?€ [0, 1[, / 6 C a if and only if it is n times continuously 
| differentiate and for all x^y ) 



Define also 

(2-27) ** = {/: (HH)7MeI ? } (a>0,.p>l). 



4 



I 



66 A. Cohen and I. Daubechies 

It is well known (and easy to check) that T +1+£ c P r n« f„ 
£ > 0. For compactly supported functions /, ^ also have ' 

( 2 - 28 ) /GC a implies far 



M» ™.H = (if) ?H . 

The infinite product (2.11) is thus divided a two parts. The first part, 

iy, since 



(2.30) 



+ 00 



[[ cos(2-y 

fc=2 



2 y 
- sin(- 



w), can 

a poiynornial expression. Indeed, since p(0) = 1 and p is a regular 
function, the infinite product generated by the second factor satisfies 

(2.31) 

Defining, for j>0, 
(2.32) 



+00 

n K2- v 
fc=i 



^ n iKrt)i- 

l<*<log(l+| w |)/lo g 2 



5j = sup 



Sp(*V 



(2.34) 



+00 

n?(rv 

fel 



i log2 ' 



<q^(i+H)/io g 2< C(1+H)i . 



and 
(2.35) 
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K«)l<0(i+M)* 



Consequently, tp is in and C a if a < N - bj- 1 for some j > 0. 
We see here that N must be large to allow high regularity since bj is 
always positive. In fact, one can prove that if the wavelet, is r times 
continuously differentiate then it has at least r + 1 vanishing moments 
(see [Mel], [Daul]), i.e. 



(%)W = 0, 



forn = 0,...,r + l and thus iV>r + l. These cancellations are 
also known as the Fix-Strang conditions [FS]; they are equivalent to 
the property that the polynomials of order JV - 1 can be expressed as 
linear combinations of the {y(x - k)} ka : However, these conditions 
are necessary but not sufficient to ensure the regularity of the scaling 
function since the effect of N may be killed by -a -large value of bj. 
Fortunately, this can be avoided b 
(and, in the biorthogonal case, additionally mo(w)). 

In the CQF-orthonormal case, a \ 
indexed by jV has been constructed in [Daul]. This construction uses 



N-l 



(2.36) %)= ^(»-l+i]^ 

(with the shorthand notation y = sin 2 (w/2)),. which is the lowest 
degree solution of the Bezout problem 

(2-37) P*)(l-# + 2/%(l-y) = l.. 
iding filters are defined by 





(2.38) 




with 






1 


(2.39) 


■ f : ■ 
; I-;. ; 





M4 = m = ?J^ 
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The Fejer-Kesz Impnta that there exists a FIR filter p w ( w ) 

JJ to (2.36 and t he conditions in Theorem 2.1 are trivially satMed 
r * " J- 1 ', f )' Por % of JV, the regularity <#) of 
the associated scaling function is approximately 0.2 Jv ffi d the exact 
asymptotic ratio between a(N) and S can he determined. Intuitively 
speaking, this means that the contribution of p N ( u ) removes "eighty 
percent of the regularity" brought by the factor (1 + e'")/^ For 
his estimation, we need to optimize the inequality (2.35), L find the 
best possible exponent for the decay of 

n.3.b. Optimal and asymptotical Fourier estimation: The role 
of fixed pomts. 



We start by defining "the critical exponent of m 0 (w)» : 
H 



(2-40) i = infi,- = iufma, _L w 
)>» j>o u-ea j log2 6 



Ten,, was proved in [Co2] that under the hypothesis W[>b(0)| = 
1 (saWed m the present case (2.39)), fl u ) cannot have a better decay 
t infim^ than If (he infimum J is attained for some finite / 
o - Oj, then this estimate is optimal. ' 

How can we estimate the critical exponent? A first method consists 
in evaluating 1, for large values of j. Weed, b is also the limit of the 
sequence bj because the boundedness of p implies bj < L + OtjlJ) 
This may however require heavy computations. 

In several cases, it is possible to use a more powerfulmethod based 
on the transformation r : u „ 2u mod-2f and the fixed points of its 
powers r», n > 0. Indeed, let Uo be afixed point of r" for n > 0 ^d 
efine,tsorbit U) = r)w „fori = 0,...,n-l. Since ft) has period 
4?r, wenave 

(2.«) ?(2 n '«i) = ?( Ui ), for all > 0 

and consequently 



(2.42) 



n 
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Letting* go to +00, this leads to 



*>-r-r kg 

nlog2 6 



for the regularity index. In fact they can do much better and provide 
optima estimate for certain types of filters. Let us consider the small- 
est or it of r different from {0}, namely the pairf 2t/3, 2r/3}. Note 
that because our filters have real coefficients, |m 0 ( W )| and \pL)\ are 



The following result 
associates the value |p(2tt/3)| and the critical exponent I 

Theorem 2,3. Suppose thtp(u) satisfies 
(2.44) 



t 2 - 44 ') WK2w)|< 



,2tt , 

VT<»<f . 



log 2' 



We already know from (2.43) that b > log |p(2ir/3)| /log 2. 
We now use the bounds on p to find an upperbound for h j > 0 We 
can regroup the factors in (2.32) by packets of one or two elements in 
order to apply either (2.44) or (2.44') on each block. Since only the last ' 
tactor can miss one of these two inequalities, we obtain 



(2.46) 
and thus, 



H 


< 


-(!) 


*=0 







J-l 



-log 



; sup [loglpj] 
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which leads to 
(2.48) 

and to (2.45). 



log 2 



The equality (2.45) means that the worst decay of occurs for 
the sequence w k = 2^/3, n > 0. This is interesting, because (2.44) 
and (2.44') turn out to be satisfied in many cases and in particular for 
the whole family of CQF defined by (2.38), (2.39). This is easy to check 
directly for small values of iV, since the inequalities can be rewritten as 

M h(y)<P^) if y < | , 



(2-49') %)W-y))<(p w (|)) if | <s/< 1 . 

The discussion for general N is more difficult and we refer'to [CC] for 
a complete proof of (2.49), (2.49'). However, a similar result can be' 
obtained in a simple way. To characterize the asymptotical behavior of 
the critical exponent when N goes to -fee, one does not need the full 
force of (2.44), (2.44'), however. It can also be derived from a weaker, 
asymptotically valid inequality, as proved by H. Volkner in [V]. 

Theorem 2.4, Let h{N) be tk critical exponent momtd to m 0 w (w) 
mi a(N) the Holier exponent of the corresponding scaling function. 
Then 



(2.50) 
and 



»-+» if 21og2 



(2.50') lim = Hm - — '±-1 = \ _ "& v „ a W t 
y ^+oo N ^ +00 N 21og2 U ' 5 ' 

PROOF. This result can be viewed as a consequence of Theorem 2.3, 
but it can also be proved directly by using some properties of P N (y) 
Let us write (2.36) in the following form: 



(2-51)- P N (i 
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From (2.36) we see that ft (1/2)- = 2 N " 1 -; since P R is an increasing 
function between 0 and 1, we have 

(2.52) P^tmaxf^qf-U^f-i. 

It is now trivial to check that (2.49) and (2.49') are satisfied if we replace 
p u{y) k j(y). The same argument used in the proof of Theorem 2.3 
leads then to 



Ml P)<giog 

but from (2.43) we get 



= 2!^ 3 



(2.54) 



' -2lb'° S 



'2JV-2W3 



JV-ll 



J-2, „ 



This proves the limit (2.50), and consequently (2.50') since the decay 
index of the Fourier transform is equivalent to the Holder exponent 
when both tend to +oo. 

The use of fixed points for optimal estimations of the spectral decay 
is thus very efficient when one is looking for arbitrarily high regularity 
since a sharp asymptotical result is obtained. For small filters, this 
method does not give a good result because the error on the exact 
regularity may have the same order as the value of the Holder exponent 
itself. For such filters, other methods, which take advantage of the 
small number of taps in the filter, can be used to derive more precise 
estimations. We now describe these methods; they are typically based 
on matrix computations. 



II.4. Regularity: Matrix based sharper estimates.. 

II.4.a. The Littlewood-Paley approach. 

We first recall some aspects^ of the Littlewood-Paley theory. Let 
7(i)-be a real-valued, symmetrical function of the Schwartz class <S(R), 
which satisfies 



(2.55) 



7(w) = 0 ifH<l/2orH>5/2, 
7(w)>0 ifl/2<|w|<5/2, 
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so that the frequency axis is covered by the dyadic dilations of 7. In- 
deed, we have 



(2.56) 0<d< £ Pu)<C 2 if b^O, 



(2.57) A;(/) = 2> 7 (2>'.)*/ <=> A ; -(/) = 7(2^.)/. 

The Littlewood-Paley theory tells us that several functional spaces can 
be characterized by examining only the I p norm of these blocks. This 
is the case in particular for the Sobolev spaces W h> and the Holder 
spaces C a , a > 0. To do this, it is necessary to change slightly the 
definition of C a when a is an integer; we shall say that a bounded, 
function / is in C n if and only if f nA belongs to the Zygmund class 
A, ie, there exists a constant C such that, for all x and t/, we have 

(2.58) \r\x±y)+r\x-y)-2r^c\y\. . 



(2.59) \Mf)\\ L «<Cr* whenj>0, 

(2.59') / is a bounded continuous function. 

Note that the choice (2.55) for 7 is arbitrary and that more general 



blocks. To derive these types of estimates on the scalingjunction y>, we 
introduce a tool which will be very useful in the bidimensional case. 

Definition 2,1. Let £ 2 [0, %\ k the space ofh-yerioiic, square inte- 
grahle functions on [0, 2?r], mi C[0, 2i] the space of h -periodic con- 
tinuous junctions. Then, for any m{w) in C[0, 2?r], we define the tran- 
sition operator T m associated to m(a>) ly 



T m : L 2 [0, 2x] — > X 2 [0, 2xJ 

:H<(H 



+ m 
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Note that when ro(u) is a trigonometric polynomial, the study of T m 

can be made in a finite dimensional space. More precisely, if we define' 

(2.61) ^a 2 )JJv-: fe.., W GC^4 

Ujv, j 

then we have clearly 

(2.62) (/, m) € [E(N h JV 2 )f implies T m f 6 E(N h N 2 ) , 

This is due to the contraction whw/2 which appears in the definition 
(2.60) of T m . If c n is the n-th Fourier coefficient of m(w), then the 
matrix of T m in the basis of the complex exponentials is given by 

W . ■ ^n = (2c 2 ,_ n ). 

The size of this matrix P in E(N h JV 2 ) isLxL with L = ft - ft + 1. 
This operator has been studied by J. P. Conze and A. Raugi and several 
ideas presented below are due to their work [CR], [Con]. We shall use it 
to derive Littlewood-Paley type of estimations for the Holder continuity 
of the scaling function. For this, we need the following result. 

Lemma 2,5, For all n > 0, 

(2.64) f(T m )7(^= f V^n^ 2 '^. 
PROOF. We prove it by induction. It is clear for n = 1 since 

/w=/;HM)HM/g4^ 

,*/2 

= 2 Kw)/(«) + m(w + + j)] <L 
Assuming (2.64) for n, we obtain at the next step, 



1 
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+ m(2~"- 1 u; + 7 r)/(2- n - 1 W + 7r)] 

/ /2 ihv 

+ m(w + 7r)/(w + 7r)](lw 

n+l 

]] m(2" i W ) /(2— ^Jdj. 



= 2 n+1 



This concludes the proof, 



We now suppose that m(w) is a positive trigonometric polynomial 
inE M = £(-M,M) and that m(0) = 1 and m(7r) = 0. Thenm 
can be factorized as 

(2.65) m(,) = cos 2W (^( w ) 

where p(w) is a trigonometric polynomial that does not vanish for w = 
7T. Note that necessarily JV < M. From this cancellation property, we 
can derive, 

Lemma 2.6. {l,l/2,...,2- 2W+1 } art eigenvalues ofT m . The tow 
vectors p j = {n%„ M ^ M , jor 0 < j < 2JV - 1 generate ft rotate 
toiici it ie/f invariant 6y T m and conteiw one eigenvector jor each of 
these 2N eigenvalues. 

Consequently, the orthogonal subspace kjinei by 

( M 



M N 

Y n\ = 0,j = 0 ? ...,2JV-.l 

n=-M J 



u ri^f invariant by T m . 



PROOF. The factorization in (2.65) is equivalent to the cancellation 
rules 

M 

(2.67) £(-l)Vc n = 0 for; =0,... ,2^-1. 

n=-M 



1 
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In particular, for j = 0, we have ■,•■■>;■ 

(2.68) Y c 2n = Y c 2«+i = » (because m(0) = 1) . 

n n 

This means that the sum of each column in the matrix of T (2.63) 
is equal to 1 and that p 0 = (1, . . . , 1) i s a left eigenvector for the 
eigenvalue 1. For 0 < j < 2N - 1 we define q, = Pj P = ( ? 7* , . . ( q ¥ )• 
we have, 3 ' J 

Thus, if I is even 

m ^E(»+|) ;<h 

and if I is odd 



that.'jj is alii 
Pi is 



! , we see 



is thus equal to 2" ; 



.Pi);=o,..,2JV-i is a 
T m and the eigenvalues are {2~ J } J=0i ... ) 2Ar_i . 

We now come back to the scaling function y>, given by the infinite 



(2.71) 



+ 00 



««) = [[m(2~* w ) 



Theorem 2.7. Let F N bt the invariant subsface of T m defined by 
(2.66). If\ is the eigenvalue ofT n restricted to F N with largest mod- 
ulus, and if |A| < 1, then, we have, with a = -log|A|/log2 (> 0), 
• ipisin C a ~ £ for a//g.>0, 
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• ip is in C a if the restriction ojT m to the invariant sub$pw F\ oj 
eigenvalue \ is purely diagonal (i.e. = 

These two estimates are optimal ify(w) does not vanish on [-7T, jr], 

PROOF. Consider the trigonometric polynomial 

(2.72) fy(«) = (l-co8w)*. 

It clearly belongs to'fjv. 
Consequently, for all n > 0, 



1/2 



(2-73) <(bf\J PMfohj 

<C{\\\U) n orC|A| n ifT m | n =AJ. 
We now use Lemma 2.5 combined with the inequality 
(2.74) C N (u)>\ when^<|«|<7r. 



J2-»f<H<J»» h»-h<\u\<2H 

y2 n 7r n 

<C C N {2-\) TTm(2-*w)d£j 
Un « £=1 

J -It 

Consequently the Littlewood-Paley blocks satisfy the inequality 
(2.75) to<Cr^ )£ >0, a = -k 



(2.75') HAjMlii < C2^', if T m | Fx is purely diagonal. 
Since ||Aj(y>)|/,<» < ||Aj(^)[ L i we obtain the announced regularity. 
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s are optimal, we u 
the inequalities which have been used. First, note that since m( W ) and 
are positive, we have IIA^)^ = |i Ay^f]^ . ■ ^ ^ ; 

Let fx be an eigenfunction in F\, E J/ A ( w )db> 0, then ' 



A(2" 



n w)]im(2-^= %fh(w)k : 
Jt=l </-ir 

If ^ /a(w) to < 0, then we replace / A by -/ A , i If jf / A ( w ) ^ = q- 



then the 



r ; see^ below after' 



?(w) does not vanish on [-x, tt], we have 

n 

(2.77) ^(w) > CjJ m(2-V) for all n > 0 and' |w | < 2 u tt . 
. fc=i 

Note that this hypothesis corresponds to the condition of Theorem 2.1 ' 
with 1 = [-*, 4 In. a more general setting, we could replace the 
integrals on [-2 n *, 2 n ?r] by integrals on Ti 
hold. Combining (2.76) and (2.77)' gives 



J~2 n Tt 



(If j ^ /a(w) = 0, then a slightly more sophisticated argument 
will do the trick. Lemma.2.5 still holds if the measure dw is replaced by 
any other measure of the type g(w) oh where g is a 27r-periodic, strictly 
positive, continuous function. We can always choose g such that 



/ A(w)j(w)4;>0; 



(2.76) then holds if fa is replaced everywhere h g{b))k. Since g is 
Vrt 1 v ^sion of (2.76) combined with (2.77), 
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Since /a has a zero of order 2JV at the origin, the function 7(1), 
defined by flw) = |/ A («)| ^(w) is convenient for the Littlewood- 
Paley analysis of Holder regularity less than 2JV. This is the case for ip 
since 2N + 1 vanishing moments would be necessary for a higher Holder 
exponent than 2iV (see [FS], [DL] or [DyL]). Consequently (2.78) tells 
us that tp cannot be more regular than C a . To prove the optimaHty of 
C a ~ £ when T m \ Fx is not purely diagonal, it suffices to replace / A by a 
function g x such that T m g x * \g\ + ji/ a with /< f 0. This leads to 

which proves the optimality of C a ~ E . 
The theorem is thus completely proved. 

Remarks. 

• The estimates (2.75) and (2.75') can be found by an equivalent 
technique, using the transition operator T ? corresponding ty the 
factor p(w) in (2.65). We simply consider the eigenvalue A p with 
largest |A| p and iterate T p on / = 1. This leads to 



/ , |w|" 

J2i-h<upit h>-h<\u\<2in 



flp(2-V 



<C(\\\ief2-W 
(orC|A p P'2^ \iyh f = \ p l) 

and thus <p € C a ~ £ with a = 2iV - log | \ p |/ log 2. This estimate is 
in fact the same as (2.75). Indeed, if p is an eigenvalue of T m in 
ffr, then its associated eigenfunction can be written as 

(2.79) /,= (rfg))\(»). 
Replacing m(w) by its factorized form in 

(2.80) ^) = /,g) m g)+/ f g + ,)m g+„) 
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we obtain, after dividing by [sin 2 (w/2) cos 2 (w/2)] N , . 

We see here that the eigenvalues of T p are exactly given by h = 
2 p. This proves the equivalence between the two techniques 
• In general m(w) is not a positive function. One can then define 
%)■= |m(w)| 2 and use the operator T M associated to M>) 
The result is an estimate of the I 2 norms of A^). Using the 
Cauchy-Schwarz inequality, we derive the following corollary, 

Corollary 2.8. Suppose that M[u) =\m(ufh s a zero of order 
Wdw = t . define \, the largest eigenvalue of T M on E N ml 
a = -kgA/(2 log2). Then, tp £ JJ«-« c M 2 ^ dere H> is 
the Soholev space of index s, The value a is attained if T M \% = 

Note that the Holder exponent has no chance of being optimal 
because we have used the Cauchy-Schwarz inequality and <p(u) is 
not a positive function. The Sobolev exponent however is optimal. 
The regularity^ compactly supported wavelets was estimated with 
this methodln [DaulL 



The transition operator plays also a crucial role in U u ,. , 
wavelet theory: we show in Appendix A how it can be used to prove 
that the families and t6Z are unconditional bases, with 
weaker assumptions than the boundedness of (1 + UW^uyi + 
|«|) imposed in Theorem 2.2. • j J ™ 

The optimal estimate for the global and local Holder regularity of 
any wavelet can be estimated by another method developed by LDau- 
bechies and J. Lagaxias in [DL]. We now recall its main points. ! 

II.4.b. The time domain approach. 

Letm(w) = Ef=o c ne mw be a trigonometric polynomial such 
that m(0) = 1 and m(t) = 0. We do not require that m( W ) be positive 
Let tp(x) be the scaling function defined by the infinite product (2.71). 
p 'It is at least a compactly supported distribution in [0, N], 
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In the time domain approach, we represent M by its "vector" 
formw(x):[0,l]-»R N 

(2.82) [w(x)] n = p(« + n-l), n = l,,..,JV. 

From the two scale difference equation (1.5) we get 

%w{2x) if x < 1/2, 
lTiw(2i-l) if i > 1/2, 

where T 0 and Tj are N x iV matrices defined by 

(2.84) (T 0 )y = cw- H Ki,j<JV, 

(2.84') = l< h j<N. 

Using the notations 

4(2;) = n th binary digit of i 6 [0,1] 

(2x if x <l/2 
TW "[2i-l if a > 1/2 (binary shift), 

we can rewrite (2.83) as a "fixed point" equation 

(2.85) «) = T<l(.)«Mx)). 



This leads to an evaluation of w(j) and its derivative by an iterative 
process. The regularity of the result depends of course on the spectral 
properties of T 0 and Tj. Note that when ro(«) has a zero of order L 
(as for the transition operator. studied in the previous section), then the 
space ft orthogonal to the vector p, = (n J ') nsl| „. )W forj = 0,...,X-1 
is invariant by T 0 and Tj. This method gives sharp estimates on the 
local regularity in x by considering the products T i:{z] • • • T w for all 
n > 0. The main result on global regularity proved in [DL; Theorem 
3.1] is the f" 



Theorem 2.9. Swppst that there exist p < 1 such that, for all binary 
sequence (dj)j^g and all m > 0, we have 



(2.86) 



\m r .T dm \ h \\<Cp m . 
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Btjine a = -log>/log2. Tkn, 

• jf a is not an integer, <p belongs to C a 

• if a is an integer, f is almost Lifschitz: for almost allxj,' 

lv— x (* + *) - V— < C|*| |log|«| | . 

Remark. 

• The "generalized spectral norm" 



(2.87) p(T h T : ) = limsup max \\T ix T di -T d \ Fr \\ l l m 



m~+oo «;' = 

j=l,...,m 



gives a sharp, estimate of the global regularity. Note that it is 
in general superior to the spectra! radius of % and T h When 
N is not too large it is. possible to compute -the exact: value: of 
p{T h T 2 ). For example, in the case of orthonormal wavelets, the 
optimal Holder exponent was found in [DL] for JV = 4, 6 and 8. 
The same evaluation becomes more difficult for larger filters. ' 
• The generalization of this approach in higher dimensions is not 
trivial. In particular, it involves nonstandard binary . expansions 
depending on the dilation matrix which is used. We describe these 
techniques in Appendix B. 



i of this review of regularity estimators, we could 
say that these three approach are complementary: the time domain 
method gives sharp results but it is only practicable for small filters, 
the Littlewood-Paley estimates can be derived for longer filters but 
they will be optimal only if m(w) is a positive function and finally, 
the Fourier approach is less precise but appropriate to asymptotical 
results on very large filters. Let us also mention that another method 
recently developed by 0. Rioul [Ri] and based on ^(Z) norms estimates 
of the iterated filters leads to interesting results; in particular, it is still 
manageable for larger filters than the time domain method of [DL]. 

' We are now ready to deal with the bidimensional wavelets. We 
start by examining the different subband coding schemes that can be 
used to build these non-separable multiscale bases. 
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III, Two channel bidimensional subband coding schemes. 



f } we shall concentrate on the dilation ma- 
trices of determinant equal to 2 or -2. In such conditions, the subband 
coding scheme that we consider split the signal in two channels (instead 
of four in the separable case) and only one wavelet is then necessary 
to characterize the detail coefficients at each scale. We first present 
a short summary of the equations satisfied by these filter. They are 
immediate generalizations of the results presented in III. 



III.1, General conditions for exact reconstruction. 

As in the one dimensional case, the scheme that we are considering 
here is based on four fundai 



• The action of two analyzing filters, one low pass 

Mo(w) = Mo(wi,w 2 ) 

and one high pass Mi(w) = M\(u h « 2 ) , 

• Decimation on each channel by keeping only the samples on the 
sublattice r = DZ 2 , 

• Insertion of zero values at the intermediate points of Z 2 /P , 

• Interpolation by two synthesis filters, one low pass 

Mo(w) = Mo(wi,w 2 ) 

and one high pass M^u) = M l (wi,w 2 ), followed by reconstruc- 
tion of the original signal by summation. 

We see here that the conditions for perfect reconstruction will not 
depend on the dilation matrix D but only on the sublattice V = DI? 
that is generated (different matrices may lead to the same T). More pre- 
cisely, there exist only two types of grid corresponding to a decimation 
of a factor 2 in Z 2 : 

• The quincunx sublattice, shown on figure 5, is generated by the 
integer combinations of (1, 1) and (1,-1). 

• The column sublattice, shown on figure 6, is generated by the in- 
teger combinations of (0, 1) and (2, 0). It is of course equivalent to 
the row sublattice, by exchange of the coordinates. 
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, The same arguments that were used in ILl.b show that perfect 
reconstruction is achieved by FIR filters, if and only if they satisfy 
(up to a shift) the following equations, which are similar to (2.13) and 
(2.14). 

t In the quincunx case, 



Mo(w)Mo(w) + M 0 (w + (7r ) f))Mo(w + (7r,?r)) = 1 



and 



, (32) m=<^^ wp+M 

4(«) = e-^ +w »)Af 0 ( W + (ir,ir)) 



Figure 5 

Quincunx decimation. 



Figure 6 

Column decimation. 
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• In the column case, 



(3.3) M 0 (w) M 0 (w) + M 0 (w + (ir, 0))M o (w + (tt, 0)) = 1 
and 



(3.4) 



&(«) = e" iwi M 0 (w + (7r,0)). 



(3.5) 



of the CQF condition (2.5). The formulas (3.1) and (3.2) become 

|M 0 («)| 2 + ^(w + Cit.t))! 3 = 1, 
Jlfi(«) = e- i ^ilfo(« + (ir,x)); 
whereas (3.3) and (3.4) become 

: i^M' + itfo («+(f,o))f = i, 

M 1 ( W ) = e- iwi M 0 ( W + (7r,0)). 

As in the one dimensional situation, we want to build from these sche- 
mes the associated scaling function which can be viewed as the limit of 
ion a 



III.2, Non-separable scaling function and wavelets. 

If c mn are the Fourier coefficients of M 0 (w), it. 
(3.7) M,( W ) = M D ( Wl , W2 ) = ^c lme - i (™»«+»»'), 



a two 



(3.8) ^) = 2^c mn ^ I -(m,n)) 

m,n 

and its Fourier transform can be expressed as an infinite product 

+00 

(3.9) fa) = [J Mo(D*"'w) 
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which is convergent, if and only if M 0 (0) = 1. 



/ if Jlf 0 (w) 

3 an FIR filter. We see from (3.9) that <j> will be highly dependent 
on the choice of D. For the same sublattice and the same filter, the 
results can be completely different for different D. The column sublat- 
tice for example is generated by both matrices D { = r 0 and 

Gi) 

Indeed, we would have 



D2 = [ 1 n ) 5 but the first one cannot lead to an I 2 scaling function. 



+00 



for all n > 0. But since ^ is compactly supported and belongs to 
£ 2 (l), it is also in L\t) and its Fourier transform should tend to zero 
at infinity. We can also remark that only the eigenvalues of D 2 have 
their modulus strictly superior to 1. " * : • 

The choice of the dilation matrix is thus very important. In fact, . 
although the equations (3.1)-(3.2) are different from (3.3)-(3.4), the 



Di such that d{!} is the column' sublattice, we can define 



(3.10) [ D 2 = PD 1 P- 1 with P : 



1 1 



Clearly, the image of f by D 2 is now the quincunx sublattice. Then, 
for any filter Mj(w) satisfying the column-CQF condition (3.6), the 
idini 



+00 



>2 is also a. seal 



(3.11) 
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! ; ' = ( j J j, we have 



WMI a +W(«+(f,f))p=W(^ I ^ + « > )P 

+ W(wi+7r,w 1 +o; 2 + 27r)| 2 = l. 

And thus I? satisfies the quincunx-CQF condition (3.5). A similar 
result holds of course if we start from two dual filters^ 1 and ^ which 
satisfy 3.3). This shows that the scaling functions associated to ft and 
are inked ^ the sim P k Nation fo(x) = ftPt). Consequently we 
can restrain our study to the quincunx case. More generally, if ft and 
ft satisfy " 1 



^2 

(3.12) 



Di = PftP" 1 



where P is a matrix having integer entries and determinant equal to 
, then we also have the same type of equivalence between the scaling 
unctions. For this reason, we shall only consider the two simplest 
dilation matrices of determinant 2, which cannot be related as in (3 12) 
since they do not have the same eigenvalues: 



(3.13) R = 



(3.13') 



1 -1 
1 1 



[Rotation of | and dilation of V5) 



1 \ ) (^Symmetry with respect to ($ + 1, 1) 
and dilation of $j , 



In both of these cases the image of f is the quincunx sublattice. The 
wavelet ^ is then defined by 

(3-14) im = M 1 (4( u ) with D = R or 5, 

where M^) is defined by (3.5) in the orthogonal case, and by (3 2) in 
the biorthogonal case where we also have a dual wavelet defined by 

(3-15) kty = Ai(u)k>) with D = R or S. 

The goal is now to design filters leading to regular scaling functions and 
wavelets. We end this section by presenting two important families of 
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filters The regularity of the associated ^, t I ^ 2 will be , , , 
it tiie one dimensional tools thai we ii 



HO. Filter design. 

III.3.a. The orthonormai case. 

Recall (see (Daulj) that in ID, the CQF filter can he designed in the 
1) For a given number JV of vanishing moments, define m, by 

(3 ' 16 > KH 2 = h 2 g)fp,^g)] 

Are P N (y) j, a polynomial, solution of the Bezout problem 
(3.17) / JMW) + (i-v)^ W = !.. 
The minimal degree choice is given by 



d solving (3,16)." 



antees that there exist 



guar- 




Unfortunately.thislastresultdoKuotgenerafiKtohigherdimen 
s«. We thus have to find other mean, to build trigonltS 
nornials wkch satisfy (3.5). One possible method i the t y £ 

S'S ; tt e ^ kft «'M^s(3,)ifand 



(3.18) 

ff oH= i(* M+Ml(w+(M) ' ¥ 'H+w,( u +(^))\ 



88 A. Cohen and I. Daubechies 

is unitary for all w. Since the product of two polyphase matrices is also 
a polyphase matrix for a third pair of filter, infinite families can be con- 
structed by multiplying elementary building blocks of the type (3.18) as 
soon as we know some simple niters which satisfy (3.5). The disadvan- 
tage of this method is that it does not furnish the vanishing moments 
in a natural way. Recall (see [Mel]) that the N times differentiability 
of the function $ implies 

(3.19) \M<C(h\ m ±k\ N +% (HhO) 

and thus M 0 (w) has necessarily a zero or order N + 1 at the frequency 
« = (tt, tt). This can also be viewed as the Fix-Strang condition (see 
[SF] ) for the regularity of the scaling function <j>. 

The simplest way to build such filters with N arbitrarily high is to 
remark that if m 0 (w) is a ID solution of the CQF equation (2.5), then 
the 2D filter defined by 

(3.20) M 0 («) = M 0 ( Wl)W2 ) = m 0 ( Wl ) 



n (3.5). It is apparently a g 
ing regular wavelets since it has the same order of cancellation in (tt, tt) 
as mo(w) in t. This allows us to build an infinite family of filters with 
s by posing 



(3.21) 



- where }mJ / (w)} N>0 is the family of filters designed in [Daul], defined 
by (2.35), (2.37) and (2.38). Note that the filter (3.21) has a unidimen- 
sional structure but since the dilation D contains either a rotation or 
a symmetry, the final analysis (using iterates of the filter) is performed 
in all the directions of the plane. In Section IV, we shall take a closer 
look at the associated wavelets and their regularity. If D = fl, then one 
can also derive another family of "almost" one-dimensional filters M 0 
from unidimensional m 0 (they get again fanned out to other directions 
by applying IT 1 ). Explicitly, 



Mp(wi,w 2 ) 



♦i 



m 0 



wi -a>2 
T 



_m o I — 2 — + 7r 



'P. 
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This construction corresponds to afflter with taps on two diagonals 

1""V A" 1 1 mi * ~" 2 + L " is *» 10 ^t 
this M 0 satisfies (3.5) if m 0 satisfies h(uf + \ n J u + t f _ j .» 

m,(0) = 1, «,(») = 0, then M 0 (i,i) = 0 Mows, so that Mi, as 
defined m 3.5) satisfies M,(0,0) = 0, as it should. One easily checks 
however, that ^( f , f ) and W)f) ^ both be - ^ 
these examples, so that the corresponding bases cannot possibly be C 1 
Only the small examples are therefore of any interest; it seems possible 
numerical experiment) to construct a continuous «S correspondirit to a 
4-tap filter in this way. 

III.3.b. The biorthogonal case. 

The filter design is clearly easier in the biorthogonal situation One 
can start from a given filter M 0 ( u ) and find the M M,(») by solving 



^ Particular we.™ look for filters which have more isotropy than 
those of the family (3.21). Here, again, tt e one dimensional theory can 
help ns to budd families of filters in a simple way. Several examples 
of real and symmetrical dual filters have been designed by the authors 
andJ.C.FeauveauinpQF], 

In these one dimensional construction the symmetry allows us to 
use the variable j = sin 2 ( u /2) and to write the transfer functions as 

M : m,(«) = p(s) and m,(w)=p(y) 
. where p and n are two polynomial satisfying 
(3.23) KsfMv) +K1 ~ 3*)P(1 - V) = 1 . 

In two dimensions, consider the variables j, = sin 2 ( Ul /2) and y 2 = 
sm (uj/2). If the filters are symmetrical with respect to the vertical 
and tbe horizontal axes, the duality condition in (3.3) can be rewritten 
as 

where 



%,l/2) ="Mo(Wi,W 2 ), P( M ) = 'M 0 
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We see that a possible choice for P and P is given by 

( 3 - 25 ) AM2) = p(aj/i+(l-a)y 2 ) 

t 3 - 25 ') HviM) =P(%+(l-a)y 2 ) 

where a is in [0,1]. For an optimal isotropy it is natural to choose 
a - 1/2; m this case the diagonals are also symmetry axes. This 
choice is known in signal processing as the McClellan transform of the 
ID filters? and p. Using the variable z - fa +y 2 )/2 we can thus write 

(3.26) M 0 (w)=p(z) an d # 0 (w) = fl*) 

where p and p are polynomials satisfying (3.24). These polynomials 
must also satisfy 

(3 ; 27) p(0) = p(0) = 1 and p(l) = p(l) = 0 

which are necessary for the construction of wavelet bases. Note that we 



(3.28) 



■!(*(?)+*(?)) 



and thus z can be regarded as the transfer function of the filter which 
computes the discrete Laplacian with the formula 

(3.29) ' A ^ m ' n = g ^ Xm,n " Xm ' l ' n " Xm+1 > n 

Since a Laplacian scheme has frequently been proposed in image 
processing to detect the edges with a maximum isotropy (see [AB], [M]), 
it seems tempting to use z or one of its powers as a high pass analyzing 
filter (and thus 1 - % as the corresponding low pass synthesis filter). 
This can be achieved in a very simple way, by a method already used 
to build biorthogonal bases in Z 2 (R). Recall that 

AM 
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is the lowest degree solution of the Bezout problem 

If we fix the reconstruction low pass as M?{u) = (1 - ( so that the 
analyzing high pass is, up to a shift, the JV-th power of the Laplacian) 
then a possible choice for the dual filter is given by 



(3.31) 



■^M = (i- 



where L is a positive integer indicating the cancellation order of M 0 at 
w - (*, 4 L has to be chosen large enough so that both functions 
m and satisfy the necessary conditions to generate a pair 



Mma* of unconditional Riesz bases (see Theorem 2 1 and 
Appendix A). We shall examine the properties of these functions and 
give an estimate of the minimal value of L in Section V. 

We have now at hand .two families of filters, orthonormal and 
biorthogonal, with an arbitrarily high number of vanishing moment, 
We still have to know if these filters allow us to build wavelet bases with 
an arbitrarily high regularity as in the one dimensional case ([Daul] 
)Co2]). As we shall see in the next two sections, the results of our' 
investigations are very surprising and show that the multidimensional 
situation contains a lot of new difficulties from this point of view. 

IV. Orthonormal bases of non-separable wavelets. 

• Let us consider the family of CQF filters defined by 



(4,1) 
with 



Mo" K«2) =. <(«i) 



and the associated sding functions for the dilations 5 and Ji , 
(«) >» = nMT(rt»), 
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(4.4) 



fc=i 



IV.l. Orthonormality of the translates. 

A first requirement is that the ^-translates of fe,s or <j> N)R are . 
orthonormal, This is a necessary and sufficient condition to generate 



Theorem 4.1. For all N > 0, the functions fc,D km orthonormal 
translates ani generate wavelet hases of the type 

v^d-h-k), jez,*ez 2 , 

roAerc d-S orR. 



Proof. By at] 

is ensured if and only if |^(w)| > Q > 0 on a c 
to [nr, tt] 2 modulo 2ttZ 2 which contains a neighbourhood of the origin. 

It is clear that M 0 w (w) vanishes only on the vertical lines wj = 
(2i + l)jr, k € Z. Consequently we see that the simple choice A' = 
[ht , tt] 2 is not convenient since for both dilations, we have 

(4.5) irW)=M) 



(4.6) ftr,,) = 0. 

Recall that in the one dimensional case, the trivial choice K = [ — jr, -ir] 
was convenient for the family m 0 N (w). Here we have to use a compact 
set K slightly different from [nr , tt] 71 so that D~>K fl {wj = (24+ l)f } 
is empty for all j > 0 and for all h in Z. This can be done very easily by 



lsoi(7r, 7r)and(-7r, nr 
them by (— 2tt, 0) and (2tt, 0) as shown in figure 7. 



One checks easily that all the sets D~>K for j > 0 are contained 
in the strip |wi| < tt - e, e > 0 where M^(w) does not vanish. 

We now have to check the regularity of the scaling functions which 
have been obtained. We shall see that the results are completely differ- 
ent depending on whether one chooses S or R as the dilation matrix. 
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Figure 7 , 

The convenient compact set K congruent to .[-tt, tt] 2 : 
1 (?r, 7r) and (-tt, hp) have been: 
wish on/{. 



IV.2. The symmetry dilation case. 



ilation matrix is S = f j ^ 



S~ l - 5 5. Since Jlf 0 w (w) = mf^), we have to consider the sequence 
{[S~H}j>o for a given u = (w h w 2 ). Clearly, it has the following 



5M«l), ^i, |( Wl + W2 ), ^i,... ) 2> 1 + W2 ),r Wl ,... 
Since T 2 = I/, the odd and the even parts are simple dyadic sequences 
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and this leads to 

(4.7) Mw) = fe(wi+w 2 )^(wi) 

or 

(is) h,s( x ) = M x i)M x i-*i) 

where y N is the one dimensional scaling function. The associated 
wavelet is defined by 

(4.9) fe>) = M 1 (w)fe )S (w) = fc(w 1 +w 2 )^(wi) 
or 

(4.10) j) N)S (x) = i) N (x 2 ) w{xi-x 2 ) . 

We see here that the scaling function and wavelet are in this case sepa- 
rable in the sense that they can be expressed directly in terms of the one 
dimensional functions <p N and fe. This separability can be explained 

by the fact that S is similar to the matrix ^ J j , which is simply a 

dilation by a factor 2 (in one) direction, followed by an exchange of the 
axes. The regularity can of course be made arbitrarily high since it is 
directly given by the Holder exponent of y N . 

REMARK. Theorem 4.1 is not necessary here to prove the orthonormal- 
ity of the translates since it is a trivial consequence of the separability 
formulas (4.7) and (4.8). 

We now consider the case of the matrix R which is by far less 
trivial. 

IV.3. The rotation dilation case. 

We now have R = J[ and IT 1 = ^ [J .The 
sequence {[iT^}^ is then, 

j(Wl+W2). \>h |(wa-Wi), -j(«l+Wj), 
-ft, ^1^2), ^l), »(«,+«,), ^ 2) ... 
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Here the first power of R~} proportional to the identity is fl" 4 = -Ij 
Consequently, it is not possible to use the one dimensional scaling func- 
tions and wavelets to express the and in a separable way. We 
first consider the case N = 1 which corresponds to the Haar filter. The 
result of the cascade algorithm with this filter shows how different the 
situation is when R is used instead of S, 



IV,3.a. The twin dragon, 

For M 0 J ( W ) = (1 + e"^)/2, the function satisfies 

( 411 )- M*)=Mk)+M&-(M)) 

and 



fc=l 

i I • By iteration of the cascade algorithm, one finds that ^ is the character- 
istic function of a well known fractal set called the "twin dragon" (see 
[K]) shown in figure 8. This set can be defined directly in the complex 
plane as 

(iI3) A = | £ "(t)" : ■ i^efo,!}"} 

and it is clear that ,6 ljfl = Xa . solves (3.41) since we have 

1 w a= (t) au (t)< a+1 ) 

f -^Aur^A + tCl)). 

I The self-similarity 'of A is thus expressed by the two scale difference 

equation (4.11), but furthermore, since the family R ( x - k)} ka > is 
orthonormal (by Theorem 4.1) and since|A| = ^(0) = 1, these inte- 
ger translates constitute a fractal tiling of the whole plane I 2 (similarly 
; i to the squares obtained in the tensor product situation with the same 
filter)., This beautiful property has been observed independently by 
: W. Madych and K. Grochenig [MG] and W. Lawton and H. ResnikofF 
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■ More generally, such tilings can be derived by considering a two 
1 5 J \c type 



(4.15) 



where D is a dilation matrix and (cj}^ are d representatives of 
2 n /DZ n (d = |detD|), This scaling function and the corresponding 
wavelet do not seem however of great interest for image processing: not 



fractal curve. Nevertheless the twin dragon is important in estimating 
the regularity (local and'global) of the wavelets with dilation matrix R, 




■0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 



Figure 8 

The "twin dragon" set A. 
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Indeed, if we want to generalize the method of [DL] (see Section II 4 6) ' 
it is necessary to consider the expansion of any point in C in terms of 
the power of (1 + i)/2 (* r 1 ), which also means that the point is 
considered as the limit of a "dragonic sequence" {A,-} >€2 with Ay c 
A H and |A ; -| = TK These "draranic expansion" technic,™ U 



Let us now examine the functions obtained with higher order filters 
which have more vanishing moments. 



IV.3.b. Higher order filters. 

We are interested in the family of scaling function <j> N)R) N>1 
Recall that in the one dimensional case, the asymptotic'result en- 
suring arbitrarily high regularity (Theorem 2.4, Section II.3.6) is based 
on the value of |m 0 (±2tt/3)| since {-h /3, 2tt/3} is a cyclic orbit of 
w » 2w modulo 2i. In the present case similar considerations for a 
fixed orbit of <* h fl w modulo hl\ lead to an opposite result: arbi- 
trarily high regularity cannot be obtained by increasing the number of 
vanishing moments. More precisely, we have 

Theorem 4.2. For oil N-> 0, the junction <j> N)R is not in C\l 2 ) , 

PROOF. This is of course true for N = 1 since we obtain the twin 
dragon. For JV > 1, we shall prove a stronger result: the decay at in- 
finity of j> N)R (w) cannot be majorated by C (which is a necessary 
condition for J> N)R to be in C 1 because it is a compactly supported 
function). For this we consider the orbit of « h Rq modulo 2t1 2 
given by the four points (2ir/5,4ir/5), (2^/5,-4^/5), (-2^/5,-4^/5) 
and (-2f /5, 4f/5). Let us denote vq = (2ar/5, 4tt/5) and v } = R\ 
One checks easily that 

( 4 - 16 ) \hM\ = C N ^ for all > 0 . 
We then have, for all iV>0, 



From the definition of m 0 w we have 



JV /27T 



(4.18) 



■©I'M-©): 



4 
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: and we know from (2.51) that 

(ii9) m<w N ~\ 4<2/<i. 

Because cos 2 (tt/5) > cos 2 (tt/4) = 1/2, we can write 



>1- 



'(31 



JV-1 



and thus, since \vj\ > 2 ; '/ 2 , 

\M v i)\ ^ ft 
>C N \vj 



i/2 



with aiv = |log[l-sin 2 (7r/5)(sin 2 (27r/5)) J /log 2. Clearly a N 
is decreasing with iV. Since on ~ 0.6115 < 1, this ends the proof. 

In fact, these wavelets do not even seem continuous although we 
have no mathematical proof for this. A simple look at the result of the 
cascade algorithm for the 4 tap filter (which corresponds to a .55 Holder 
continuous one dimensional wavelet) shows how chaotic the functions 
h,N can be (figure 9). The design of FIR filters leading to regular 
wavelet bases with R as the dilation matrix seems to be a difficult 
problem. Using a polyphase component approach M. Vetterli and J. 
Kovacevic ([KV], p. 32) have constructed a filter for which the result 
of the cascade looks continuous but no infinite family with arbitrarily 
high regularity has been designed so far. 

The main difficulty which makes this design impracticable is the 
absence of the Riesz lemma in more than one dimension and thus the 
impossibility to start by designing the square modulus of M 0 (w) in an 
appropriate way. Apart from this problem, the CQF filters (in particu- 
lar the family (3.21) that we have introduced) cannot be symmetrical. 
We must keep in mind that one of the interests of the quincunx grid 
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decimation is to have a more isotropic analysis; this is c 
the filter coefficients are themselves symmetrical a 
vertical and diagonal directions. 

These two reasons encourage us to construct biorthogonal b 
wavelets from dual filters for which the Riesz lemma is not necessary 
and linear phase can be achieved. 



1 bases of 




Figure 9 

ion of the scaling function fan, 
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V.Bi 



Let us recall the family of dual filters introduced in III.3.b. It is 
based on the variable 



(5.1) 



where L is still to be fixed. 

A first remark is that the action of the dilation matrices R and 
S on the variable z are equivalent. This is due to the fact, that z, is 
invariant if we exchange wj and w 2 or if we change the sign of one of 
■ these variables. We shall thus consider a dilation matrix D which can 
be equal to R or 5. To express its action on z we still need the two 
variables 

(5-3) «,-*(&) and S2 = sm J (|) . 



z = ^(!/i + !/2)^2 = ^i + t/ 2 -2 M2 ) 

D 1 



M « s ,(*li(l-ld+*ltf-«0) = S<l!+»i) 



We shall start by studying the scaling function f associated to the 
filter Mj(w) = l-j, because it is the elementary building block for the 
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V.l. The quincunx Laplacian scheme. 

The coefficients of Mj(w) are centered around the origin and have- 
the following form: 



(5.4) 



/ 1 \ 
14 1 

\ 1 / 



Note that this is the simplest symmetrical filter (with respect to the 
horizontal, vertical and diagonal directions) which satisfies the cancel- 
ation MjfTr, f ) = 0. To estimate the decay of ^(w) 



J] cos(2" fc w) 



sinw 



used in the one dimensional case. Note that (5.5) is based on the iter- 
ation of sinw = 2.sin(w/2) cos(w/2). Unfortunately, similar relations 
do not exist in the bidimensional case for the dilation matrix D. In 
particular the infinite product 



&(«) = J] Mo(D- J W ) 



has no simple expression and one checks easily that, unlike (5.5), it 
does not have uniform decay at infinity. Indeed, let us consider the sets 
{(27T/5, 4i/5)} and {(2*/3, 2ir/3), (2jt/3, 0)}. These are two cyclic 
orbits of w h Dw modulo Iff and modulo the exchange of coordi- 
nates and sign changes which do not affect the variable z. Consequently, 
if we define Vj = D j (2k /5, 4* /5) and ^ = j0i(2ir/3, 2x/3), we 




(5.7) k(vj)«C 



2 (?r/5) + cos 2 (27T/5) 
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with 

(5.9) a v = ~ 
log2 



; 2 (7r/5) + cos 2 (2tt/5) 



2.83 



2 ; 13 



2 2.68^ a, . 



Still we would like to find a global exponent for the decay of ^(«) 
at infinity. For this we shall introduce an "artificial" function which 
will play the same role as cosw in (5.5). We define 

(5.11) CH = - Vy U = 



Contranly to %), q w ) i s not a trigonometric polynomial, but it 
is a bounded regular function which vanishes at the point (», tt) with 
the same order of cancellation as Mj(w). Moreover, it satisfies by 
construction 



+ 00 



(5.12) JJ^P 



.js 2 [sin 2 (m/2) + sin 2 fo/: 



The decay of this infinite product is now uniform and, for this r„ 
son, C(u) will play an important role in the construction of our dual 
bases. F 



Proposition 5.1. The decay ojfa(u) at infinity is controlled by 

M lliHI<c(i+H)- 2 . 

Furthermore, this exponent is globally optimal, i.e. there exists a se- 
quence {(4 >0 such that lim^ +00 ^1 = +oo and \^)\ ~ C\^\~ 2 . 

yi = an 2 (wi/2) and y 2 = sin 2 (w 2 /2) we 



can rewrite C(w) as 



(5.14) q w ) = yi+^- 2 viSfa = (i ~ yi)y 2 +(i - y 2 )y! 
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We thus have 



1/1+1/2 2 

_ (1 - yi)(y 2 ~ y 2 )(yi — y 2 ) • 
%+») . 

_(yi-!/2) 2 



>0. 



Thus Mi(w) < C(u) and by (5.12) < C(l + H)~ 2 . To prove 
that this exponent is optimal we consider a small vector p j. 0 in R 2 
and define 

(5.15) «j^(f,f) + p, 



(5-16) h(»t) = J[Ml(&-%*) + D- k P ) 
Let us divide this product in three parts 

Wwi)'=[i((x,T) + H)] 
"J'-I 



n^(^(',')+ir*p 




•K(M+n)] 

V«)C(j). 

One checks easily that ^(tt, x) / 0 and thus, for j large enough or 
sufficiently small />, we have 0 < C, < < 1. It is also clear that for 
1 < H j - 1, Ml ,)) = 1 and that for m, Jlf }(% T )+ 
ff) > 1 - C \\a\\ for a small enough, with C > 0. Consequently, if p has 
been chosen small enough, 1 > B(j) > f]^ [1 - C2"< ||p||] > C 2 > 0. 
Finally since (tt, tt) is a second order zero of M 0 (w), the third factor 
satisfies 

m ™ 3 \\ef = c 3 \\D-j ? t< C (j) 

<c 4 rvi 2 = 2-^|rf. 
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This shows that behaves like V * |wj|" 2 when j goes to -j-co 
and the proposition is proved. 

Note that from the decay of ^(«) we cannot even conclude that 
it belongs to L\l 2 ) or that fa(x) is a continuous function. Yet both 
are true; we are going to prove this by the Littlewood-Paley method 
explained in IL4.a. The filter Mj(w) and the scaling function ^(w) 
are particularly well adapted for this approach since they are positive 
sojhat the regularity estimation is optimal (because ||Aj(^)||ioo ~ 
\\£j(4>i)y\ see Section II.4.a). 

Proposition 5.2. The optimal glohal Holier exponent for fa(x) is 

PROOF. We consider the transition operator defined by 

(5.19) TF(Du) = Ml(w)F{w) + M 0 ] ( w + (x, i))F{w + (x, x )) . 

As in the one dimensional case T can be studied in a finite dimen- 
sional space but this subspace cannot be defined as simply as £?( ^ N 2 ) 
in (2.61). One way of finding an invariant subspace is to apply T to the 
constant 1 and then iterate it on the characters e i (* lUl+ * lW *) which are 
obtained until a stable set is attained. With Mq corresponding to (5.4), 
this subspace is trivial, since T\ = 1. Lemma 2.5 then guarantees the 
integrability of j> h hence the continuity of <j> h To estimate the Holder 
exponent of ^ we need a larger subspace, which we obtain by iterating 
T on 1 and on cosw! + cosw 2 . The size of the matrix representing the 
action of T on this subspace can be seriously reduced by exploiting the 
symmetries, i.e. the invariance under u : <— > -w h w 2 <— * -u> 2 and 
V\ *— >w 2 . 

Using the subspace E generated by the basis 

(5.20) ci = 1, e 2 = coswj +cosw 2 , e 3 = cos^ +w 2 ) + cos(w 2 -wj) 
we obtain the following matrix 



(5.21) 



/l 1/2 0\ 
0 1/2 1 
\0 1/4 0/ 
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which has the eigenvalues {l, (1 + V5)/4, (1 - The two last 
eigenvalues correspond to the subspace ft C E defined by 



(5.22) 



ft = {F(u)tE\ F(0) = 0} . 



Similarly to the one dimensional case, we iterate T on the positive 
function e } - \e 2 which is clearly in ft and this leads us to 



where Aj^i) is the Littlewood-Paley block corresponding to the re- 
gion ^([-tt, x] 2 )/IT 1 ([-x, tt] 2 ), situated at a distance 2>/ 2 of the 
origin. Consequently, if we define. . 



(5.24) 



2 , I- A 



then it follows from (5.23) that 

(5.25) {n\u\f^)a\l 2 ) and ^(i)€C°(R 2 ), 

Consequently ^ is Holder continuous with regularity 0.61. ■■ 

This property appears in the graph of ^ on figure 10 (obtained 
by the cascade algorithm) which presents a smooth aspect with several 
pointwise cusps. ^ Note. that this regularity is not: sufficient to derive a 
better decay of ^(w) than |w|" 0,61 ; Propositions 5.1 and 5.2 are thus 



Remarks'. 
• Note that, since we have 



26) X>) + M!(« + M) = i, . 

we can derive the L l convergence of the truncated products 
= IIj=i M)(0~Mxd»([-m] 2 )(w) with the same method as in the 
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£^ + 2k) = l 



ier and Dubuc [DDI It " ****** 
™ c \ uu \- « explains the four 



wps surrounding the center 

(-1,0) which are visible on fW e m Z" ^> 



, iW = 0aret 




Figure 10 



( c i>e2,V4,e s ) > 



C08 ^ w ( w i+« J ),coe(u 1 - Wj )) 

/ 2 A t-\ 0 0 \ 

» 1-A A 0 2 

0 A 2 0 

OA 0 0 0 
\° 0 1-A o o/ 



W J ty SS attaiaedforA = 0orl Note W ) ™' mde « f V 
convolution product »M - » i } COrres P t>nds '» the 
introduced in lV Q ** *" A Wn **» 



the function ffl ,1 7p" l ~ ^ W w «ffl 
2/1 + 2/2 



2/i > 2/2, and 

2/1+2/2 "t 1 -^) 
' 2/i + 2/2 

^ 2/2 > 2/1 .On the other, hand 



^0 



(5.30) 
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to majorate for 2»> < |,| < frDIt we on j y need to v 
rate the j first factors in (5.30). Since R rotates by tt/4, half of the 
factors can be majorated by C(u) and the others by 2C(w) This 
leads to ; 



|^)|<C^H)/lo g 2 jj C(r , w) 

KK21og(l+| w |)/lo g 2 



and thus 



iM<C(l + |w|) 



It is easy to check (in a similar way as for £(«)) that this esti- 
mate is optimal. An immediate consequence is that the Fourier 
transform of the twin dragon characteristic function X a satisfies 



(5.33) 



XaM<C(1 + H)-V2 



which was not obvious since we did not have a formula similar to 
(5.5) for fo. 

We now return to the construction of our biorthogonal bases and 
attack the problem of obtaining isotropic wavelet bases with arbitrarily 
high regularity. 



V.2. Biorthogonal wavelet bases with arbitrarily high regu- 
larity. 



We now consider the whole family of filters 



L>0 



i by (5.1) and (5.2). 
A first remark is that the regularity of the functions increases 
linearly with N. More precisely, since 

(5.34) 



feW = [*)%(*) 



we can use the characterization of the optimal decay exponent for &(w) 
established in Proposition 5.1 to estimate the regularity index «{N) of 
fe(a;). This leads to ' 

P-35) 2N-2<a(N)<2N 



Nob-separable m mmm waveiet eases 



and thus to 
(5.36) 



. Jim _ n 
»-4» JV ■ ' 



The estimate (5.35) is of course more interesting for large values of JV 
For N = 1, we have seen that a ~: 0,61, 

the svmmetr.es reduces the size of the matrix to 9 x 9, Analysing th 

rr'/Tf tbt ^ is inC "* a , 2.93. The functb 
h - <h * h looks very smooth indeed on figure 13. ■ 

For JV > 3 the matrix becomes too large to tackle by hand. In all 
cases the regukrty of the wavelet ^ will of course be the same as 
ha of The problem is now to find the appropriate dual function 
for the analysis. More precisely we want to design the filter 



osing the nmnber I in such way that the hypotheses of Theo- 
rem 2. ,n its hdimensonal generalization) are satisfied, i.e. that we ' 



) <C(l + (t 



h-l-t 



£>0. 



Tc i show that such a choice is possible for any value of JV («.' for M 
arbi ranly regular synthesis function), we need an asymptotical 1 result 
of the same nature as Theorem 2.4. We want to be sure that the 
regu taMg action of the factor (1 - ,f m compemate the ^ 
enect 01 r Nj . L it i ]S large enough. 

Using a similar .approach, we consider the simplest fixed point of 

l\f*n* h f MmM "^^^feexcLgeof. 
.Tksfixdpomtis^ = (2,/5,4,/5) which corresponds 

W 2 0 = Z[Uq) = 5/8. 

We now decompose M 0 w,i into three factors, by introducing the 
function C(u) defined by (5.11): 6 

with 
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We already know from Section IL3.b that 



(5.40) 



PN{z)<{izf A ih> 



from the Bezout relation (3.30), we also have 



(5.41) 



1-z 



Consequently, we can roughly majorate Pjv(z) by 



(5.42) ftr(*)< 



min{- — , max{4z,2}} 



if z € [0, 1] . 



Defining JTfw) = min{l/(l-z), max{42,2}}andG(w) = H{w)Q(u), 
(5.39) leads us to 

(5.43) ^)<[C(4\(K4m7 • 

We are now facing a similar situation as in Theorem 2.4 where we had 
shown that the function g(y) = max{2,4y} = h(u) satisfied 



(5.44) 



f.W=*(v)<*(3/4) if 2/ < 3/4, 

[^)W = ff(y)j(4|((l-y))<b(3/4)] 2 if 3/4 < y < 1 . 



In the present case, although we do not dispose of any simple mathe- 
matical proof, numerical evidence shows that we have 



(5.45) 



(5.46) 



|G(w)G(Dw)<[GW J or if not,. 



lH(u)H(Du)<[Efaf or if not, 



{H{u)H{Du)H(D 2 u)<[H(u<>f 
These two statements are illustrated respectively in figures 11 and 12. 
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(a) 




o h 

a) Graph of G : {y h y 2 ) = mu{G{w)G(D^ [G( Wo )f } - [<Z( Wo )f 
(b) 




o h 

b) Graph of G 2 (y h .y 2 ) = mu{G{u)G{Du)G{D 2 u\ [G(« 0 )] 3 } - [G{u a )f, 
(c) 




c) Compared supports of Gi(y h y 2 ) and G 2 {y h y 2 ) 
Figure 11 
Graphic proof of (5.45) 
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a) Graph of y 2 ) = max{%)# (Dw), [%)] 2 } - [%)] 2 
(b) 

2.0 \ 




o h 

b) Graph of ftfa, y 2 ) = max{tf( W )jp w )tf (i) 2 w ), [%)]»} - [%)f 
(c) , 




c) Compared supports of fl^, y 2 ) and ff 2 (y 1 , t/ 2 ) 
Figure 12 
Graphic proof of (5.46) 



! : 



r 
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On a) and b) of each of these figures we have plotted the functions 
and 

for F = G and H (the coordinates are {y h y 2 ) g [0, if). On c) the 
supports of a) and b) are shown to be disjoint regions in [0, l] 2 . 
We now estimate ^(w). From(5.39) and (5.43) we get 

+00 

fc>)=n[^«)] L %(^) 

Jc=l 

<^(i+Hr 2i n w - 

l<K21og(l+| w |)/lo g 2 

<^(HHr 2L [ n 

l<K21og(l+| w |)/log2- 
TT l ] N 

n *p-m • 

KK21og(l+>|)/log2 



) to divide these products in groups of two of 
three factors which satisfy one of the inequalities, this leads to 

(5.47) l N>L (w) < C(l + |w|)" 2 W 1 * G M + iviogif( Wo ))/iog2 , 



(5.47') fe, L (o;)<C(l + H) 2i M+2^ 

with 

Fortunately a < 1. This means fe^i) can be made arbitrarily regular 
by choosing I large enough. In particular, (5.38) will be satisfied if we 
have 

(5.48) %-l) + 2W < -1. 
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The smallest I such that U» and fi^ generate unconditional mul- 
tiscale bases is therefore given asymptotically by 

(5.49) L(N) * JL = lQ g5 - lo g 2 

1-a logl6-lo g 15 " 14,2jV ' 
This asymptotical estimate is moreover optimal. Indeed define w • = 
Vo* = ^(2ir/5,4f/5). Because of the fixed point property of «* 
. we clearly have ' 

M fV^c[M„%)]'', Chr7 
with 



(5.51) 



log2 



Rom the definition (5.2) of M^ 1 ', we get 



>C 



■ 3 a . . »« 



8j 2 



-err 



and thus 



7 <C + 2X»-2JVM? 

log2 l og 2 

= C + 2L(l- a )-2^. 

It follows that the estimate (5.49), if true, is certainly optimal. While 
we expect (5.45), (5.46), hence (5.49), to be true, we have unfortunately 
no rigorous proof. However, we can prove inequalities which are slightly 
less strong than (5.45), (5.46), leading to a non-optimal but rigorous 
estimate for L(N). More precisely, we can prove that ft = U d 2 
can be split up as ft = ft, uft 2 Uft 3 , with 

|GW<£ w€fl,,' 
M (?(a;)G(D w )<(2 wGQ2) 

lG( w )G(^)G(D a «)<f» W 6ft 3) ' 
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with (12- . 9588 <1, resulting in (5.47') with 

hi 

log2 

If we use the crude estimate H{w)< 4 for all w 6 K t] 2 , correspond- 
ing to ^ = 2, then this leads to 



. L(N)~JLn« 32.959 JV; 

this factor is about twice as large as in (5.49), The 
this estimate is in Appendix C. •■= iU - ■■ * •■ 

All these results can be summarized in the following-theorem. 

Theorem 5,3, The family of iul filters {M^\ M^ NL> , 
generates horthogonal hases of compctlfsu^oHd mvtltis with ar\i- 
tmily high regularity. For large values ofN,the Holier exponent of 
fa(x) is equivalent to 2N ani the minimal choice for L is asymptoti- 
cally proportional to N , 

P-53) . L(N) ~a , 

iwiU4.215'<JC< 32.959. 

Here the upper bound on I is not tight, and we expect £ =. 14,215 
to hold, as indicated: above. 

REMARK. By taking L larger than L{N)^ can also be made arbi- 
trarily regular. However, in many applications such. as coding, approx- 
imation, data storage and compression, we do not really care about 'the 
regularity of the analyzing functions j> and only the synthesis func- 
tion i> and <j> have toibe smooth since this property is important for the 
cascade-reconstruction algorithm. This justifies the choice of the min- 
imal value L(N) such that the families {2>( 2 j) N)L (D>x -.k)}^ kg? 
and {Vl 2 i NtL (Dh - k)} mk£P are unconditional dual bases of 
L (I). Recall that the existence of frame bounds is essential for the 
stability of the subband coding scheme. 

We end this section by taking a closer look at the size of these dual 
filters. 
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V.3. Size and optimal implementation of the dual filters. 

The asymptotical ratio L(N)/N a 14.2 is big in the sense that the 
filter M 0 JVi(w) may have a very large number of taps. More precisely, a 
polynomial P(z) of degree p corresponds to a filter with p 2 + (p + 1) 2 
nonzero coefficients. For example, if N = 3, 

is according to (5.49) a polynomial of degree p = 2L( jV) ~ 87 in z. 
Consequently it is the transfer function of a filter with approximately 
1350 taps! 

It seems thus that the dual filter is, even for small values of N, 
much too large for a realistic implementation. This is not quite true for 
several reasons. 

First, one can factorize the polynomial P^l{n)( z ) ^ express 
M,f ,L as a product of p monomials in z, By applying successively these 
monomial filters instead of using directly their product, the number 
of multiplications per sample in the filtering process is reduced from 
order p 2 to p. Note that this complexity reduction associated with 
the factorization is due to the multidimensional situation and does not 
occur in the ID case. 

Second, the filter corresponding to the variable *, it. the laplacian 
discrete scheme, has coefficients c 0|0 = 1/2 and c lj0 = c_j 0 = c 0 1 = 
. c 0 -i = -1/8. It can thus be implemented by using binary shifts instead 
of multiplications. This is very important since a binary shift is usually 
performed 10 times faster than an addition and 100 times faster than a 
multiplication in most processors. This shows that only the additions 
count here. If t is the time for one multiplication, each monomial filter 
will generate one sample in approximately 3t/5 and the same operation 
will take 3pf/5 for the whole filter. For iV = 3 and p = 87, this 
corresponds to the complexity of a 52 tap filter which is much more 
reasonable than the first estimation. 

Finally, for small values of N } it is clear that the asymptotical esti- 
mate (5.49) of L(N) is far from sharp, just as, in ID, the asymptotical 
estimate on regularity of Section II.3.b was ill-suited to small filters. A 
better estimate for L(N) can be found by checking that the optimal 
decay exponent for #(w) is exactly determined by the value of M^' L at 
w 0 = (2x/5, k/5). More precisely recall that we have 

<- L ( u ) = [CHf%H. 
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For the small values of N and L considered below, one can check by 
the same graphical arguments that the inequalities (5.45) or (5.46) are 
also satisfied by % i.e. 

(5 54) (%( W )%(^[%MP or if not, 
Wiv,^)%(D W )5 NjL (D 2 W )<p iV)L ( Wo )] 3 . 

In order for (5.38) to be satisfied, we therefore only need 
(5.55) M 0 %)<| , 

and this will be sufficient for these small values of N and L for which 
(5.52) holds. Using the definition (5.2) of M^ L we obtain 

• forJV = l,L(l) = 3, 

• for JV = 2,1(2) = 12,' 

• forJV = 3, £(3) = 22. 




Figure 13 

. The scaling function fa (= ^ * fa ). 
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Clearly these estimates are much better than (5.49). Finally, L{H) can 
be even more reduced, for small values of tf, if an even sharper criterion 
that the frequency decay (5.38) is used to ensure the existence of frame 
bounds. We show indeed in Appendix A that the spectral analysis of 
the transition operators T 0 and % corresponding to the functions |M 0 1 2 
and |M 0 | 2 can be used to derive both the frame property and the L 2 
convergence needed to have a pair of dual Riesz bases. In [CD] we prove 
that this criterion is sharp so that the value of L(N) is here optimal. 
Unfortunately the matrices of T 0 and % can be very big, even for small 
JVandl. 

For N = 1, we now obtain L(JV) = 2 so that the two filters M\ 
and M 0 12 are of small size. We show on figure 14 and 15 the scaling 
functions and wavelets obtained from such a choice. 



Appendix A: A sharp criterion for frame bounds. 

We want to give here a better result than Theorem ft to char- 
acterize the dual filter pair (m 0 , m 0 ) which lead to biorthogonal Riesz 
bases of wavelets. The method that we show here uses the transition 
operators associated to the positive functions |m 0 | 2 and |m 0 | 2 (see Sec- 
tion II.4.a). 

First, recall that the <p, (p, $ and i are defined by 



(A.1) 



#(w) = JJ m 0 (2 _ *w), i(2u) = mi(w)0(w), 

+00 



As mentioned in Theorem 2.2 the duality relations (2.19), (2.20) and 
the decomposition formula (2.21) are ensured as soon as the partial 



Ww)=n^(2- fc w)x h2 « jr)2njr ]( w ; 

n 

Vm(«)=n *«( 2 "* W )Xl-2»*,2-ij(w; 
Jb=l 

converge in L 2 [l) respectively to p(w) and p(w), 
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The main difficulty is then to obtain the frame bounds A, B A 
and B all strictly positive such that for all / in X 2 (R), 



(A.2) 



m\ 2 <i\(fM<B\\ft 
^i/ii ! <E'i(wi)f<%i ! 




It is sufficient to obtain the two upper bounds of (A.2) because 
the lower bounds are then obtained by using (2.21) and the Schwarz 
inequality which give 

/ yi/2/ y/ 2 

M !/f< \Z\ifM\- 

U* / U / 

In [CDF] we used the following assumption 

( Ai ) rtw)p + ^ u )p.<C(l.+ Hr 1 ^ 

which can also be formulated with <p and (p instead of ^ and i Here 
we shall prove the ^convergence of { K ^ n>0 ^ the frame in- 
equalities (A.2) using weaker assumptions. More precisely let T 0 and 
T 0 be the two transition operators associated to the functions |m 0 | 2 
and K| ,.as defined in Section II.4.a/ They both operate in two 
spaces of trigonometric polynomials E N and fy. We have proved in 
Lemma 2.2 that the subspaces F N = {/ 6 E N , : /(0). = 0} and 
^JV - {/ ■€ Efj : /(0) = 0} are invariant under the action of T 0 and 
T 0 . The following result gives a sharp characterization of the dual filter 
pairs associated to biorthogonal wavelet bases. 



Theorem A,l. Let \ (mytctivtly A) he tk largest eigenvalue of T, 
(respectively., T 0 ) in tk suh m F N {m^ktly, F$< Then if \\\ 
mi \\\ are hoth strictly inferior to \, tk functions j> and j defined hy 
(A.l) generate korthogonal Riesz km of wavelets t^ m , 

PROOF. We shall prove here that this condition on the eigenvalues of 
To and T 0 is sufficient to obtain biorthogonal wavelet bases. In fact it 
is also a necessary condition. This result is detailed in [GDI; 
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We first show that ip and j are in L% As in Theorem 2.7, we 
apply T 0 n to the function d («) = 1 - cob « which is in F N and by using 
Lemma 2.5, we obtain 

A~!«<M<s"t h'-',<M<iH 
f' 

<C C,(2-" u )|v„( w )| 2 (t, 

J-2H 

because ( 7 + 1)/2 > 7 , Since we also have ( 7 + 1)/2 < 1, it follows that 
the dyadic blocks in the Littlewood-Paley decomposition of ip satisfy 
the inequality 

( A - 5 ) \Mv>)k<C2-<> forsomeOO. 



1 are even better than L 2 : They belong to a 
Besov space (c L\l)) for some e > 0. We shall use tl 



to prove the frame inequalities. Similarly (p and $ belong to B^ for 
some e > 0. To prove the L 2 convergence of the sequence ip n to we 
remark that since m 0 (0) = 1, for a in ]0, *] small enough we have 

(A.6) | w |<a implies |$w)|>C>0. 



( A -7) ») = n m o(2"Mx[-2-a )N M.' 



It is clear that f n (u) converges pointwise to $w), but (A.6) also im- 
plies |#;(«)| < |fl w )|/C for all n > 0. By the Lebesgue dominated 
convergence theorem we get 

(A.8) limK-* J = 0 . 

We now use the hypothesis on the eigenvalues to evaluate the L 2 norm 



J J\u\> a 



f.; 
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■ 1 A 
n-+oo , 

Consequently % converges to f in V(l) wi the same y ds for . 
and (p. Yn 

upper frame inequalities in (A.2). We 
lemma. 1 

Lemma A.2. Let j k a jmim in L 2 (l) such tki for some 0 > 0, 

M EW w + 2 ^)l 2 " a <Ci 

Then, for all f mL 2 (l\ v ; 

( A - n ) E i(wi)i 2 <ci^ii/i 2 . 

Let us first assume that this result is true to conclude the proof of the 
theorem We thus have to check that there exist a O 0 such that 
(A.9) and (A.10) are satisfied for i> and i . 

To check (A.10), we define I; = -2%] U [2\ 2^], 



For j < 1, v 



[w).at the origin to obtain 



(A.12) ™\M<CV for ; < 1 , 



We i^ 2 Ce kn ° W ^ ±2>7r) = ° since * h) = 0for ^ z W 



max 
we/; 



<J[ 
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1/2 



The first factor can be majorated by 2 _i > because we have proved that 
^bngs to £■ . The second factor is finite since it is proportional 
to j W(x)| dx and ^ is a compactly supported £ 2 function. Conse- 
quently 

(A.13) max |^)| <C2^ 2 f or j> 2 . 

Combining (A.12) and (A.13) we see that (A.10) holds for all O 0 
since we have 



i<l j>2 /J 



We now 



some a > 0. 
= 0 for all k G Z, we can derive 



fib 



fa 



We already saw that the second factor was finite (in the proof of (A.10)). 
The first factor is also finite for o small enough. Indeed, using \j> 6 
and the Holder inequality, we obtain 



Jl i Jlj 



(2 i+1 f) ff 



■ T 
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<C2^- 2j M) ;;i 

ti iU5 nave to cnoose (j > 0 such thaf /r 9c/i n « . 

Shoe the s^ 2 Iwd&f, 0 '"; << 
the theorem is proved modulo Lei A 2 1 .^™* t * 



= S:j 0 r, EA ri («+2ftr)j^+'2ftt) 



1/0 >/a- • 



•1^ + ^)1^/2 \ 



2_i / 2r Iv^ 
%]' (.g l/( 2 "* J '(^ +^ 7 r))|2 |^ (w+ 2 



W«)P«ti. 



Sunumg pn all the scales; eZ^and using (A.10), we grt 
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Appendix B. Dragonic expansions. 

In this Appendix we want to show how the one-dimensional tech- 
niques in [DL] can be extended to multidimensional situations. As an 
example we discuss the two-dimensional case, with the dilation matrix 

A first multidimensional extension of [DL] can be found in [Mo], 
Even though he looks at general matrices, Mongeau effectively reduces 
his analysis to pure dilations by considering the smallest n such that 
D = D n is a multiple of the identity, and rewriting (by iteration) 
the two-scale equation so that it involves only D. This procedure can 
drastically increase the number of different terms in the equation. We 
choose here to work directly with D = R itself. 

. When the two-scale equation is one-dimensional, and the dilation 
factor is 2, the regularity at x of the function <j> solving the equation is 
regulated by the binary expansion of x (for dilation factor p, the same 
role is played by the p-ary expansion). Moreover, R and in particular 
supp^ is tiled with integer translates of the interval [0,1], which can 
be viewed as the set of numbers equal to the decimal part only of their 
dyadic expansion; if N such tiles are needed to cover the support of <j>, 
then the two-scale functional equation can be rewritten as an equation 
for an iV-dimensional vector-valued function involving two matrices T 0 
and 2\. The spectral properties of T 0 , 2\ then determine the regularity 
of ij), both local and global [DL]. 

In the two-dimensional case with dilation matrix jR, the role of 
elementary tile is now played by the twin dragon set A. It is defined 



(B.1) 



ha 2 -. 



P;€^ = Z 2 /iS 2 = {(0,0) 5 (l,0 



Under the standard identification of E 2 with C, with (x, y) ~ x + iy, A 
can also be written as 



(b.2) a = L e c: z = gd^i±iy where(i . = Oorl 



NON- SEPARABLE BIDIMENSIONAL WAVELET BASES 127 

This set A is compact; has fractal boundary, is selfsimik, and its Z J - 
ranslates hie the plane. The indicator function of A is the solution to 
tne two-scale equation 

4>(x) = <j>(Rx) + #Jk-(l,0)) 

(see [GM]). A is called the twin dragon set [K]. We shall give the name 
dragonic txpmmm to expansions of x or z as in (B.1) (B 2) Note 
that (as m the binary case) some points may have two different dragonic 
expansions, e.g. .01000... = ((1 + i)/ 2 f = f /2 = .1011111... (This 
example also illustrates that addition follows rules very different from 
the binary case, since .0100 . . . + .0100 . . . = .1111 . „ ) 

Suppose we are interested in various regularity properties of L 1 - 

1T1 Ann A ■ 



(B.3) 



where A is a finite subset of l\ Such solutions are. uniquely defined 
up to normalization and have necessarily compact support. One can 
determine the minimal set r C l l so that R'\T + A - L) C r- then 
mf C U^(A U). The equation (B.3) for * can be rewritten by 
dehning the |r|-dimensional vector v(x) by 



we have 



where d^x) is. the first digit in the dragonic expansion of x, and rx is 
the point obtained by dropping ^(x) from the same dragonic expansion 



pi 

Equation (B.4) can be recast as 

(R5 ' »(*) = T iM ;v(rx t 

where (% = c aHj ffl 
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d a setup analogous to that of [DLJ. The question 
is now whether the proof techniques of [DL] still work in this case. The 
answer is basically yes. For instance, we still have 

Theorem B.l. Assume that the c k in (B.3) satisfy 

n n 

Then t x = (1, is a common lefteigenmtor o/T 0 , T x vftji 
e^enWve 1 /or both matrices, Dejine E l to be the one-iimensional 
subspace orthogonal to e h If there exist \<\,C>Uso that 

( B - 6 ) \h-TM\<c\ m 

for all possible d, = 0 or 1, all m € N, then the L l -solution <j> to (B.3) 
is EoUer continuous with exponent a = | log A|/ log 

This is the analog of Theorem 2.3 in [DL]. Two different strategies of 
proof are given in [DLJ. The first one involves piecewise linear spline' 
approximants; this technique would be hard to generalize here because 
of the fractal boundary of our domain building blocks A + h. A second 
strategy, which does not use splines at all, but leads to longer proofs, 
is explained in the Appendix in [DL]; this strategy generalizes to the 
present case. The main point we have to check to make sure the proof 
carries over is whether elements that are close necessarily have dragonic 
expansions with the same starting digits. In the one-dimensional, bi- 
nary case, if two dyadic rationals i, y are closer than 2~ m , \x-y\ < 2" m , 
then x and y have binary expansions with coinciding first m digits. (If 
e.g. x < y < x + 2~ m , then the expansion "from above" of x - ending 
in all zeros - has the same first m digits as the expansion "from below" 
of y - ending in all ones.) This is crucial in the proof, and allows to 
extract Holder continuity from the condition (B.6). We therefore have 
to check whether a similar property holds in the "dragonic" case. 

^ By analogy we shall call irapnk rationals all the points in A for 
which a terminating dragonic expansion can be written. Typically drag- 
onic rationals also have other, non-terminating dragonic expansions. 
For each dragonic rational x the terminating expansion is unique; we 
denote its digits by (fj(j)j' 6 N. 

Let us also introduce the notations fl 0 , R h 

% = Ry, R iy = Ry + (1,0), 



1 
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ot % = fy + d(l,0), mty=0orl. ; .. 

Take now a fixed dragonic rational i, and assume that M - 0 
for i > J. All the j 6 A that have the same first ] digit, & ; < r 
constitute a little dragon Aj(i) themselves, " ; ' ' 

itself is the image of {0, 0) under the same map R]\ '. '■ ■ Br 1 The 
■ set A is^iled by 2^ little dragons of the S ame size Jl,, all Wlates 
of A; For every such little dragon, we call the point corresponding 
to (0,0) the_ bottom and the point corresponding to (0, 1) (the only 
other pent m Z^fl A) the "top". H x is a dragonic rationai withi 
most N nonzero digrts, then i is the bottom of AjM for all j > JV 
(But note that the, "orientation" of A,(»), as indicated by the line 
connecting bottom and top, changes with ;!). It Mows that , is on 
the order of these A;(i), If * is not at the edge of A itself, then there 
must exrst another little dragon A,(,) so that , is the top of A M 
7f;.^ofaJlthe2'possible dragons A,). Since theto 
1 of A is given by the expansion .111111..., we can therefore fin ' 

^f mc W fo ".«%malones,andwith.the S a« 
Arst digits as y, 

4W = <?M fori<;, dj(x) = i {ori>J: 

We have seen how to obtain the two expansions for a dragonic 
rational * We now want to show that if : another dragonic rational y 
■ close to i then at least one of its expansions starts with the same 
digits as one of the expansions for t. Define 

e = max{r:B((0,0); r)cAU(A-(0,l))), 

where %;A) is the open Euclidean ball centered at , with radius \ 
buppose l is a dragonic rational with <fc) = 0 for j > /. Take m > J 
and consider the set ' 

4 = (j6A: |}-i|<,)2-"/ 2 }. 

abilities: either i is on the ba 
en 
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has no common interior points with A, so that B m c A m (x), and the 
terminating expansions of allyefl B have the same first m digits 

set is then a little dragon A m ( 2 ) of which x is the top. In this case 
B m C A m (a;) U A m (z), so that every point y 6 B m has a dragonic 
expansion with either the same first m digits as <?(*) (if y g A m (x)) 
or as <?(«) (if y 6 ^(z)). This is the main ingredient needed to make 
the proof of Theorem 2.3, as sketched in the Appendix in [DLL work in 



«mic 



: ™v**iw» ror x Qoes not lead to inconsistencies 

for the definition of If d »(,) = 0 for j > J, then d\x\ dHx) are 
linked by v ' 



^E d iW^M = ^4wruo) + r>,i) 

i=i i=i 

for > J arbitrary. One can then compute v(x) in two ways,' using' 
the two expansions. The following computation shows that they lead 
to the same result: for h £ T 



/jv 



E c KH^( a )(i,o)-i 1 %- 1+ d§(x)(i 1 o)-i 2 • ■ • c %+W(i,oW, 



E C WK*)(l,0)-m 1 • •• c Hm w . 1 +<i|i(i)(l 1 0)-m JV 

E C «l(*)(l,0)-mr^^^ 
mi...mjv ' 

The reader can now check that the proof in [DL] indeed carries 
over to prove Theorem B.l. Similarly, one can prove differentiabil- 
ity of I under stronger conditions on I 0 , T h similar to Theorem 3.1 
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in [DL]. Finally, ^the same techniques can also be used for local reg- 
ularity estimates, but these are a bit more tricky, and require fur- 
ther study of the properties of dragonic expansions. In practice, the 
matrices %\ Ev Tfa are often too large to permit a rigorous esti- 
mate of A in (B4 However, \ is bounded below by the quantities 
Pftii "' r imk) 1/m , and this leads to upper bounds for the Holder 
exponent a. 

Examples. 

1. g{x) = h(& + (l,0)) + s (&) + 1^.(^0)) 
The solution to this equation is the convolution xa * Xa, where 
Xa is the indicator function of the dragon set A (see also the 
second remark following Proposition 5.2). In this case T has 10 
elements. The largest spectral radius of % is obtained for 
d = 0, />(T 0 |jSi ) = .847810 . . . , corresponding to a lower bound 
A > p{%\ El ) in (B.6), or a Holder exponent a < .47637 . . . Via 



derives that this value is a lower bound. This global H61der expo- 
nent is attained in dragonic rationals, in particular in (0,0). 
Note that when M 0 is positive, as in this case, the transition oper- 



that the matrix representing T is in fact a submatrix of T 0 , so that 
it is not surprising that they have a common eigenvalue! 

2. ^)r^«&)+^&-(l,0))+M(ito-(-l J l))+fc3«&- 
(0,1)), with hg = .506970418225, h : = -.207072424345, h 2 = 
.493029581775, k = 1.20707242435. This is an example from the 
family described at the very end of Section III.3.a 
thonormal wavelet basis. In this case \T\ = 14; L_ 
been chosen so that p(%\ El ) a pp^) ~ .714. 1 
mations to i seem to suggest that j might be continuous, but we 
have no proof. If it is, then its Holder exponent is bounded above 
k hWikf^lbifiz log(.90649)/ log V5 - .28327 . 

| Appendix C. Proof of the inequalities (5.52) for 

The function G is defined as 
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2W1+W2 jm-wj 



2 



with 



h(f) = M 

I 4 ' 1/2 < * < 1 . 



G(f «), where %, U2 ) , either («, +U2 , «, or ( U] U[ ^ 
«!). (Since G is invariant for the interchange of «,,«,, it does not 
matter which definition of fl is taken, D = or D = 5.) To prove 
these inequalities it is convenient to use different variables, 

We then have 



Moreover, 

W = 2M, p(D w ) = 4(6 2 -p). 
As « ranges over [-*, f ], (,, p ) fill out the domain A defined by 
A = {(s,p): 0< 5 <1, max{0,2 5 -l}<p< 5 2 } . 



' new vanables , ^ therefore want to study n(s p) 
iM and ^,p) ^(, jp)) ^ M)j for all ( ^ 

wnere U is defined by 

^M = M = (2(^-P), 4( 5 2 -p)). 



Note that D maps A twice onto itself (both A n {s < 1/2} and A n 
{s > 1/2} get mapped to all of A). Moreover D has~one fixed point, 
(«o,Po) = (5/8, 5/16), corresponding to ^.po) = 15/16. 
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We shall prove that A = AiUA 2 UA 3 , where 
(CI) 'iy)<t on A!, 

(C.2) i*MB(w))<-<* onA 2 , 

(c.3) onA 3. 

Thev 

«1. 



Ai - ! »)€A:; (1-^3 if «<l/2, 



P<S»- 6 2 (1~ 5 ) if S>l/2 



Since 



(C.4) 



= if 5 <l/ 2) ^ if 3> 1/2 , 

" : 2($-p) 3-p 



j|(a,p) <aon Ai 



we automatically nave 



By the definition of ij and 5, we have to distinguish four different 
regions when studying j) 2 (s,p) = ??($,p)?](l)(s,p)): 




2(1- p) 4(s-2s 2 +p) 

233(1-5)' 

s-p 

I S 2 (l-s)^ S\l-S) 

s-p 5 — 2s 2 + p 

'(1 - s)(l - J) _ 8/(1 - s)(s - p)(l - - p)) 
s-p s+p-2s 2 

if s > 1/2, p< s-1/4 



if s < 1/2, p< «-l/4, 
if s < 1/2, p>s-l/4, 
if 3 > 1/2, p>s-l/4, 
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We define A 2 by 

A 2 = j( 3) p)eA: S <l,p>Jl-.I 



Up)eA: 5>»,p>l,8 s -.81 
= A 2) i U A 2)2 . 
Since s(\ - s) < 1/4 for all s G [0, 1], we have 



2(5 -p) ~ 4( S -2 S 2 +P ) 
on all of A 2| i. Since moreover p > (1 - l/2a) s, we have 



»fe(w)< 



4 1 



2a 



<<* 2 



on Ao 



On A 2|2 fl {(s,p) G A : p > s - 1/4}, one easily checks that 
s\l-s) 



s-2s 2 +p 



satisfies d f rj 2 f 0 everywhere. It follows that rj 2 achieves its maximum 
on the boundary of this domain, given by the three pieces p = s\ with 
l/2<K.9,p = s - 1/4 with 1/2 < s <. 7, and p = 2cWwith 
.7 < 5 < .9 . One easily checks that the maximum value of ^ on this 
boundary is .9. 

Similarly one checks that i/ ? achieves its maximum on A 2 _ 2 D {(s,p) G 
A : p < s - 1/4} on the boundary of this set; again this leads to 

!»<.». 

It follows that 



,(C.5) 



%(«,p)<.9 = a 2 on all of A 2 



It remains to determine an upper bound on ?| 3 (s,p) = ?](s,p)??(2)(s,p)) 
# 2 (s,p))onA\(AiUA 2 ) = {( s ,p) : 2 S -Kp< s 2 ,p>. 
5 - 2s 2 )(l - «)/a, p < 1.85 - .81} Since s - 2s 2 (l - 5 )/<* is strictly 
increasing, we have A \ (Aj U A 2 ) C A 3 = {(s,p) : 2s - 1 < p < 
5 2 ,Pi = 1.85i - .81 < p < 1.85 - .81}, where 4 is 1 
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tos-2* 2 (l-s)/a = 1.85 - .81. In A 3 one 



\i = A 3 n{p > Pl , p> 2 5 -l, p < 2( s ~l/4) 2 }, A 3)2 = A 3 'n{p> 
2(5 -l/4) 2 ,p< 5-1/4}, A 3 , 3 = A 3 n{p>s-l/4,p,>2(s-l/4) 2 "} 
and A 3 , 4 = A 3 fl {p > s - 1/4, p < s(s - 1/4) 2 }. On A 3)1 , A 3)3 and 
A M one checks explicitly that fyi/ 3 ^ 0. On A 3)2 , the exact expression 
for rji is too complicated, but one can replace it by an upper bound, 



h(m>) 



= j^- gH 2f(l-j) 
s-p s-p l-p 

Js 2 {\-s) s 2§ 2 (1 - 1) 
s-p s-p J-p . 

_ 166 2 (1 -6)1( 1 -g) , _ 

- ~ = 

s-p 



i / 0 on A 3)2 . It follows that .7/3 on 
A 3 is bounded by the maximum of ^ on the boundaries; of A3 1, A3 3, 
A 3)4 and of t/ 3 on the boundary of A 3)2 . Explicitly, for all (s,p),G A 3 , 

(C.6) <%('i>Pi) = •88145650226... 

This numerical upper bound is larger .than (.9) 3 / 2 ; it follows therefore 
from (C.4) and (C.5) that we have proved (C.1)-(C.2) for ■ : 

. ■! ( = tt,Pi)] 1/3 = .958812370442... '■; 
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